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Abstract 

We  present  a  time-dependent  semiclassical  transport  model  for  mixed-state  scattering  with  thin 
quantum  barriers.  The  idea  is  to  use  a  multiscale  approach  as  a  means  of  connecting  regions  for 
which  a  classical  description  of  the  system  dynamics  is  valid  across  regions  for  which  the  classical 
description  fails,  such  as  when  the  gradient  of  the  potential  is  undefined.  We  do  this  by  first  solving 
a  stationary  Schrodinger  equation  in  the  quantum  region  to  obtain  the  scattering  coefficients.  These 
coefficients  allow  us  to  build  an  interface  condition  to  the  particle  flux  that  bridges  the  quantum 
region,  connecting  the  two  classical  regions.  Away  from  the  barrier,  the  problem  may  be  solved  by 
traditional  numerical  methods.  Therefore,  the  overall  numerical  cost  is  roughly  the  same  as  solving 
a  classical  barrier. 

We  construct  numerical  methods  based  on  this  semiclassical  approach  and  validate  the  model 
using  various  numerical  examples.  In  the  one-dimensional  case,  we  use  a  finite- volume  method  that 
extends  the  Hamiltonian-preserving  scheme  introduced  by  Jin  and  Wen  for  a  classical  barrier.  In 
the  two-dimensional  case,  we  consider  a  mesh-free  particle  method  that  can  be  computed  efficiently 
and  that  may  be  extended  to  higher-dimensions.  The  semiclassical  transport  model  is  verified 
numerically  by  examining  the  convergence  of  the  Schrodinger  and  the  von  Neumann  equations  to 
the  semiclassical  limit  for  several  examples.  Finally,  we  examine  an  extension  of  the  model  to 
coherent  dynamics  necessary  for  periodic  crystalline  and  mesoscopic  scale  quantum  barriers. 
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Chapter  1 
Introduction 

In  this  work,  we  develop  a  semiclassical  model  of  particle  dynamics  in  the  presence  of  a  thin 
quantum  barrier.  Thin  quantum  barriers  include  quantum  dot  and  quantum  wire  structures,  res¬ 
onant  tunneling  diodes,  thin  films,  and  interfaces  between  two  dissimilar  materials  such  as  p-n 
junctions.  Simulation  of  electron  distributions  in  the  presence  of  such  barriers  is  important  to 
the  understanding  of  the  behavior  of  plasmas,  semiconductors,  and  modern  electronic  devices.  Ad¬ 
vances  in  nanoscale  materials  fabrication  technology  have  prompted  the  need  for  efficient  numerical 
simulation  of  quantum  structures.  However,  simulation  is  difficult  when  the  system  reacts  over  dif¬ 
ferent  length  and  time  scales  since  the  smaller  scale  usually  drives  the  accuracy  and  consistency  of 
the  solution.  Even  when  only  interested  in  the  macroscopic  behavior,  one  may  be  forced  to  resolve 
the  microscopic  dynamics.  Correspondence  principles  allow  us  to  extract  macroscopic  behavior 
from  microscopic  dynamics  in  terms  of  a  weak  limit.  When  the  scales  act  over  several  orders  of 
magnitude,  the  numerical  solution  to  the  problem  at  the  smallest  scale  becomes  computationally 
intractable.  In  these  cases,  one  often  relies  on  a  multiscale  approach  to  provide  a  numerically 
efficient  solution. 

An  example  is  the  modeling  of  electron  transport  in  nanostructures,  such  as  resonant  tunneling 
diodes,  superlattices  or  quantum  dots,  where  quantum  phenomena  in  localized  regions  of  the  devices 
cannot  be  ignored.  While  one  can  use  quantum  mechanics  in  the  entire  region,  it  is  clearly  more 
computationally  efficient  to  take  a  multiscale  approach  using  classical  mechanics  in  the  rest  of  the 
device  by  using  a  domain  decomposition  technique.  Such  a  model  was  introduced  by  Ben  Abdallah, 
Garnba  and  Degond  [6,  7,  8].  In  this  model,  interface  conditions  connecting  the  classical  and  the 
quantum  regions  were  used  to  couple  two  classical  regions  with  a  quantum  region. 
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This  work  is  an  extension  of  the  Hamiltonian-preserving  finite- volume  method  introduced  by  Jin 
and  Wen  [21,  22]  for  solving  the  multi-dimensional  classical  Liouville  equation  with  a  discontinuous 
(but  classical)  potential.  The  idea  there  was  to  build  an  interface  condition,  such  as  used  in  [7], 
that  properly  incorporates  transmission  and  reflection  information  at  the  barrier  into  the  numerical 
flux.  This  produces  a  scheme  that  connects  momenta  (velocities)  on  both  sides  of  the  barrier  via 
the  Hamiltonian  preservation  principle.  Such  a  method  is  stable  in  both  l1-  and  Z°°-norms  under 
a  hyperbolic  stability  condition  and  captures  sharply  the  weak  semiclassical  limit  of  the  linear 
Schrodinger  equation  or  geometrical  optics  through  the  barrier  or  interface. 

The  quantum  barrier  that  separates  the  two  classical  regions  differs  from  a  classical  barrier 
in  that  a  quantum  wave  can  tunnel  through  a  barrier,  be  partially  transmitted  and  reflected  by  a 
barrier,  and  resonate  inside  a  barrier.  The  method  proposed  in  this  thesis  is  to  solve  the  Schrodinger 
equation  (either  exactly  or  numerically)  inside  the  quantum  barrier  in  order  to  generate  transmission 
and  reflection  coefficients,  and  then  use  that  information  in  the  interface  condition  to  solve  the 
classical  Liouville  equation  through  the  barrier,  in  the  spirit  of  the  Hamiltonian-preserving  method 
of  Jin  and  Wen.  When  the  quantum  barrier  is  thin  (on  the  order  of  a  de  Broglie  wavelength),  solving 
the  time-independent  Schrodinger  equation  suffices.  Thus,  the  first  step  is  merely  preprocessing. 
Once  the  transmission  and  reflection  coefficients  are  generated,  the  time  marching  is  based  on 
classical  mechanics.  Hence,  this  approach,  which  efficiently  handles  a  thin  quantum  barrier,  has  a 
computational  cost  similar  to  a  classical  simulation  in  the  entire  device. 

The  primary  focus  of  this  dissertation  is  the  development  and  numerical  implementation  of  the 
semiclassical  model.  While  a  model  may  be  rigorously  defined  mathematically,  we  are  often  unable 
understand  its  intricacies  without  computer  experiments  and  simulation.  Mathematician  Peter  Lax 
once  stated  “It  is  impossible  to  exaggerate  the  extent  to  which  modern  applied  mathematics  has 
been  shaped  and  fueled  by  the  general  availability  of  fast  computers  with  large  memories.  Their 
impact  on  mathematics,  both  applied  and  pure,  is  comparable  to  the  role  of  the  telescopes  in  astron¬ 
omy  and  microscopes  in  biology.”  [27]  Therefore,  a  secondary  focus  of  this  thesis  is  the  simulation 
of  the  semiclassical  Liouville  equation  and  the  von  Neumann  equation  for  several  examples.  In  this 
manner,  we  are  able  to  both  validate  and  verify  the  model  and  also  illustrate  particle  dynamics  in 
a  variety  of  environments. 
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In  Chapter  2  we  introduce  the  semiclassical  model.  The  chapter  provides  a  brief  review  of  the 
underlying  physics  and  derives  the  semiclassical  limit  of  the  Schrodinger  equation.  The  correspon¬ 
dence  between  classical  and  quantum  mechanics  is  discussed.  The  formal  model  is  developed  in 
terms  of  the  weak  Hamiltonian  property  and  an  interface  condition.  Finally,  we  examine  limitations 
of  the  model  in  terms  of  entropy  and  time-irreversibility. 

In  Chapter  3  we  propose  the  implementation  of  the  semiclassical  model  and  its  numerical  dis¬ 
cretization  in  one-dimension.  The  one-dimensional  interface  condition  is  derived  using  the  transfer 
matrix  method  and  the  model  is  employed  using  a  finite-volume  method.  We  present  four  numeri¬ 
cal  examples  to  verify  the  numerical  method  and  validate  the  semiclassical  model.  Our  numerical 
results  indicate  that  the  model  correctly  captures  the  solution  of  the  Schrodinger  equation  in  the 
entire  domain  in  the  limit  of  the  vanishing  scaled  Planck  constant. 

In  Chapter  4  we  extend  of  the  numerical  implementation  of  the  semiclassical  model  to  two 
dimensions.  The  two-dimensional  interface  condition  is  derived  using  the  quantum  transmitting 
boundary  method  and  the  model  is  employed  using  a  mesh-free  particle  method.  We  present  two 
numerical  examples  to  verify  and  validate  the  semiclassical  model  and  the  numerical  method.  The 
results  indicate  that  the  model  correctly  captures  the  Schrodinger  solution  in  the  limit  semiclassical 
limit. 

In  Chapter  5  we  discuss  corrections  to  the  model  for  thin  barriers  to  extend  it  to  mesoscopic 
and  periodic  barriers.  This  chapter  attempts  to  circumvent  the  shortcomings  of  the  model  and 
provide  a  direction  for  future  research. 
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Chapter  2 

Semiclassical  Model 


2.1  Correspondence  between  classical  and  quantum  mechanics 

2.1.1  From  classical  to  quantum  mechanics 


A  typical  problem  under  consideration  is  particle  flow  in  a  plasma  or  through  a  solid-state  device 

over  a  macroscopic  scale.  If  the  potential  is  sufficiently  smooth  we  may  describe  non-interacting 

particle  dynamics  in  phase  space  classically  as  the  Hamiltonian  system 

dx  V  ,  „  dp  ,  , 

-  =  -  =  VpH(x,p),  -£  =  -VxV  =  -VxH(x,p)  2.1 

at  m  at 

where  x{t)  E  is  the  particle  position,  p(t)  E  is  the  momentum,  m  is  the  effective  mass 
and  V(x)  is  a  time-independent  potential.  The  Hamiltonian  function  H(x,p)  represents  the  total 
energy  of  the  system 


H(x,p)  =  f-  +  V(x)=E. 
2m 


(2.2) 


One  may  introduce  a  probability  distribution  of  particles  f(x,p,t )  in  phase  space.  By  requiring 
that  the  probability  be  conserved  along  the  particle  trajectories  (the  Liouville  condition),  one  has 


With  the  help  of  equation  (2.1),  one  gets  the  classical  Liouville  equation 

=  {H,  /}  =  Vp/  •  VXH  -  Vxf  •  VPH 


(2.3) 


where  {  • ,  •  }  is  the  Poisson  bracket.  Alternatively, 

jrJ  +  -  •  Vx/  -  VxV(x)  •  Vp/  =  0.  (2.4) 

at  m 

By  considering  the  zeroth-order  moment  of  f(x,p,t),  one  obtains  the  probability  position  density 
in  physical  space 

P(x,t)  =  /  f(x,p,t )  dp , 
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which  serves  as  a  primary  observable  for  the  comparison  of  the  model. 

When  the  potential  fluctuates  rapidly  over  a  short  distance  or  the  particles  impinge  on  a  sharp 
jump  in  potential,  the  classical  description  fails  to  capture  the  quantum  wave-like  nature  of  the 
particle  and  the  Liouville  description  produces  an  incorrect  solution.  In  particular,  the  classical  Li- 
ouville  equation  does  not  model  barrier  tunneling,  probabilistic  partial  reflection  and  transmission, 
or  resonance  which  are  crucial  to  the  behavior  of  many  modern  electronic  devices. 

By  considering  Dirac  quantization,  one  has  the  formal  correspondence  between  the  classical 
quantities  and  the  quantum  operators 

x  — >  x,  p  — >  -ihV,  and  E  — >  (2-5) 


where  h  is  Planck’s  constant.  Using  this  quantization,  one  obtains  the  Schrodinger  equation  from 
the  classical  Hamiltonian  (2.2) 

ih ^  =  H^=  (2.6) 

which  describes  the  time  evolution  of  the  probability  amplitude  ip(x,t]  x,p )  initially  centered  at  x 
with  an  initial  energy  state  E  =  E[(x,p).  The  square  of  the  magnitude  of  the  probability  amplitude 
p(x,t )  =  \il>(x,t)\2  gives  the  position  density  in  physical  space. 

Instead  of  considering  a  pure  state  system,  one  may  also  consider  a  mixed  state  system  for  which 
the  initial  state  H{x,p)  of  the  particle  is  given  in  terms  of  a  macroscopic  statistical  distribution 
f(x,p).  Define  the  density  matrix  as 


p(x,x',t)  =  /  f(x,p)'i/j(x,t-,x,p)'ip(x',t;x,p)dxdp. 

J«.d  JRd 


(2.7) 


The  time  evolution  of  the  density  matrix  is  found  by  taking  the  partial  derivative  of  equation  (2.7) 
with  respect  to  t.  By  using  the  Schrodinger  equation  (2.6)  and  the  hermicity  of  Hamiltonian 
operator  H,  one  obtains  the  von  Neumann  equation 

d  (  ft 2  \ 

ih—p(x,x',t)=  f  —  [Ax  -  Ax>]  +  V(x)  -  V(x')  J  p(x,x',t).  (2.8) 

The  formal  correspondence  between  the  Liouville  equation  and  the  von  Neumann  equation  is 
seen  by  replacing  the  Poisson  bracket  in  equation  (2.3)  with  the  commutator 


{HJ} 


ih  1[H,P \ 


6 


giving 

ih^p=[H,p\=Hp-pH.  (2.9) 

which  is  equivalent  to  equation  (2.8). 

The  von  Neumann  representation  may  be  thought  of  as  the  fundamental  description  of  quan¬ 
tum  mechanics  [9].  By  taking  f(x,p )  =  8(x  —  xo )5(fi  —  po)  in  (2.7),  the  density  matrix  reduces  to 
p(x,x',t)  =  il>(x,  t;  xo,po),il>(x' ,  t\  xo^po)  and  the  physical  observables  of  the  mixed  state  von  Neu¬ 
mann  equation  correspond  to  those  of  the  pure  state  Schrodinger  equation.  In  this  manner,  the 
Schrodinger  equation  is  simply  a  limiting  case  of  the  von  Neumann  equation.  By  taking  the  diagonal 
of  the  density  matrix,  one  gets  the  position  density  in  physical  space 


p(x,x,t) 


f(x,p)\ip(x,  t ;  x,p) | 2  dx  dp. 


2.1.2  Semiclassical  limit:  quantum  to  classical 

Consider  a  characteristic  length  and  time  scale  L8x  and  L5t  where  Sx  is  the  natural  length  scale 
such  as  a  de  Broglie  wavelength  5x  =  h/p  for  some  momentum  p.  By  rescaling  x,  x'  and  t 

xt-^-x/LSx,  x'^-x'/LSx,  tt-^-t/LSt 
in  the  von  Neumann  equation  we  have 

is-^p(x,x',t)  =  ^--^—[Ax  -  Ax/]  +  V(x)  -  V(x')\  p(x,xr)  (2.10) 

where  the  dimensionless  scaled  Planck  constant  e  =  [mL(8x)2  /  8t]~xh  and  the  effective  mass  m 
has  been  nondimensionalized.  Solving  the  Schrodinger  and  von  Neumann  equations  numerically 
presents  several  difficulties.  The  de  Broglie  wavelength  must  be  resolved  numerically  to  ensure 
correct  physical  observables  of  the  solution.  Typically,  this  requires  that  the  mesh  size  Ax  =  0(e) 
or  even  o(e)  with  a  similar  constraint  on  the  time  discretization  At  [4,  33].  When  £  is  small, 
computation  is  expensive  since  we  need  to  use  0(Nd+1 )  operations  to  compute  the  Schrodinger 
solution  and  0(N2d+l )  operations  to  compute  the  von  Neumann  solution  where  N  =  0(e_1)  is  the 
number  of  grid  points  in  each  space  dimension.  Because  of  such  reasons,  semiclassical  methods  are 
important  for  the  solutions  when  £  <C  1. 

A  typical  path  to  the  derivation  of  semiclassical  limit  is  through  the  WKB  approximation. 
However,  the  WKB  approximation  to  the  Schrodinger  equation  fails  to  capture  multiphase  infor¬ 
mation  beyond  caustics  [19,  40].  An  alternative  method  is  to  use  the  Wigner  transform,  the  Fourier 
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transform  of  the  density  matrix, 

W(x,p,t)  =  j77-Td  [  P(x  +  \ey,x  -  \ey,t)eTlv'y  dy.  (2.11) 

\^)  J»d 

By  applying  the  transform  to  the  von  Neumann  equation  one  has  the  Wigner  equation  [43] 

r\ 

--w  +  —  •  vxw  -  e£w  =  o 

ot.  m 

where  the  nonlocal  term 

0£W (x,  p,  t )  =  — l—j  [  -  [V(x  +  \ey)  -  V{x  -  \ey)\  W{x,y ,  t)e~tp  y  dy 

jRd  £ 

with 

W(x,y,t)=  I  W(x,p,t)eip'y  dp 

Jmd 

being  the  inverse  Fourier  transform  of  W(x,p,t).  The  expression  for  QeW  may  also  be  expressed 
using  the  Wigner-Moyal  expansion 

°°  / _ i\n  l  e  \ 

Q?w  =  •  VPW  +  V  7  vln+1v  •  Vln+1W. 

x  P  (2n  +  l)!  x  P 

71—  1  X  7 

When  the  potential  V (x)  is  sufficiently  smooth,  one  recovers  the  classical  Liouville  equation  in  the 
limit  as  £  — >  0  [12,  31] 

%  +  --yxf  ~VXV  -Vpf  =  0.  (2.12) 

ot  m 

However,  the  classical  limit  is  not  valid  at  the  discontinuities  of  the  potential  [3,  36,  38],  where 
the  potential  behaves  as  a  quantum  scatterer.  In  the  case  of  a  quantum  barrier,  we  may  consider 
a  multiscale  domain  decomposition  approach  for  a  solution  [7].  In  the  next  section,  we  present  a 
semiclassical  model  of  a  thin  quantum  barrier  with  the  mixed-state  dynamics. 

2.2  Particle  behavior  at  a  quantum  barrier 

To  model  quantum  dynamics,  we  consider  a  top-down  multiscale  approach  by  considering  the 
quantum  effects  as  local  corrections  to  the  global  classical  particle  dynamics.  In  order  to  isolate 
and  simplify  the  problem,  we  make  the  following  assumptions/limitations: 

1.  The  effective  width  of  a  barrier  is  0(e).  On  the  classical  scale,  this  means  that  we  may 
approximate  the  barrier  as  having  zero  width;  on  the  quantum  scale,  this  means  that  we  may 
typify  it  as  a  single  scattering  center  and  we  may  neglect  particle  dwell  time  in  the  quantum 
region  in  the  semiclassical  limit. 


2.  The  distance  between  neighboring  barriers  is  0(1)  and  hence  each  barrier  may  be  considered 
independently. 

3.  The  change  in  the  potential  W(s)  is  0(1)  except  at  quantum  barriers. 

4.  The  coherence  time  is  sufficiently  short  and  therefore  we  may  neglect  interference  away  from 
the  barrier. 

Naturally,  one  would  like  to  be  able  to  treat  a  wider  class  of  problems  including  periodic  crystalline 
domains  and  mesoscopic  barriers  for  which  e  is  nonvanishing.  We  will  examine  corrections  and 
extensions  to  these  simplifications  in  Chapter  5. 

We  begin  with  the  Hamiltonian  system  discussed  in  Section  2.1 

^x  =  VpH(x,p),  =  -VxH(x,p). 

Let  a  bicharacteristic  of  the  function  H(x,p)  be  the  integral  curve  <p(t)  =  (x(t),p(t)).  Note  that 
ip(t)  may  not  be  defined  for  all  time  f  G  R.  When  H{<p(t ))  is  differentiable, 

jtH(v{t))  =  jtx-VxH  +  jtp-\7pH  =  0  (2.13) 

from  which  it  follows  that  the  Hamiltonian  is  constant  along  any  bicharacteristic  <p(t),  i.e., 

H(tp(t))  =  const.  (2-14) 

Condition  (2.13)  may  be  interpreted  as  the  strong  form  of  the  conservation  of  energy,  while  con¬ 
dition  (2.14)  may  be  interpreted  as  the  weak  form.  If  the  potential  V(x)  is  discontinuous  or  not 
defined  in  some  region  Q  €  Md,  the  Liouville  equation  fails  to  have  a  global  solution  since  VXV(Q ) 
is  undefined. 

The  key  idea  behind  Hamiltonian  preserving  schemes  [21,  22]  is  to  (a)  solve  the  Liouville  equation 
locally;  (b)  use  the  weak  form  of  the  conservation  of  energy  to  connect  the  local  solutions  together; 
and  (c)  incorporate  a  physically  relevant  interface  condition  to  choose  the  correct  solution.  Let  £ 
be  the  locally  defined  set  of  bicharacteristics  of  the  function  H(x,p).  By  requiring  the  Hamiltonian 
to  be  constant  along  trajectories,  we  create  an  equivalence  class  of  bicharacteristics  [ip\  ={</?*€: 
£  |  H(<p*)  =  H{p)  }. 

Generating  a  global  bicharacteristic  is  a  matter  of  connecting  equivalent  bicharacteristics  at  the 
barriers.  Consider  the  incident  and  scattered  trajectory  limits  (x{t~),p(t~))  and  (x(t+),p(t+))  on  a 
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barrier  in  one-dimensional  physical  space — the  two-dimensional  case  will  be  considered  in  Chapter  4. 
From  equation  (2.2)  the  scattered  momenta  are 

p(t+)  =  — p(t~ )  (2.15a) 


for  reflection  and 


p{t+)=p{t  )y/l  +  2 m[V(x(t~))  —  V{x{t+))]/\p{t~l)\2  (2.15b) 

for  transmission.  Unless  |p(t_)|2  <  2 m[V (x(t+)  —  V(x(t~))],  for  which  the  transmitted  momentum 
is  imaginary,  the  conservation  of  energy  does  not  tell  us  which  of  these  two  bicharacteristics  a 
particle  should  physically  follow.  In  order  to  resolve  the  nonuniqueness,  we  require  an  additional 
interface  condition  which  we  derive  from  the  Schrodinger  solution  across  the  interface.  By  inter¬ 
preting  a  wave  function  as  a  statistical  ensemble  of  a  large  number  of  particles  [35],  we  have  the 
interface  condition 

f(x{t+),p{t+))  =  R(pR(t~))f(x(t+),pR(t~))  +  T(pT(t~))f(x(t~),pT(t~))  (2.16) 

where  T(p)  denotes  the  probability  of  an  incident  particle  being  transmitted  across  some  region, 
R(p)  denotes  the  probability  of  an  incident  particle  being  reflected,  and  the  incident  momenta 

PR.{t~ )  =  - p(t+ )  and 

Pr{t~ )  =  p{t+)V1  +  2 m[V(x(t+)  -  V(x(t-))]/\p(t+)\2 . 

come  from  equations  (2.15b)  and  (2.15a).  The  multiple-dimensional  interface  condition  is  an  ex¬ 
tension  of  the  one-dimensional  interface  condition  (2.16)  and  will  be  discussed  in  Section  4.2. 

By  considering  the  time-reversibility  of  the  scattering  process,  we  may  formulate  an  alternative 
but  equally  valid  interface  condition 

f(x(t~),p(t~))  =  R(pR(t+))f(x(t~),pR(t+))  +  T(pT{t+))f(x{t+),pT(t+))  (2.17) 

where  the  scattered  momenta  pR{t+ )  ans  PT{t+)  are  functions  of  the  incident  momenta 
PR.{t+)  =  ~ p(t~ )  and 


PT{t+)=p(t  )v/l  +  2 m[V(x(t~)  -  V(x(t+))\/\p{t~)\2 
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To  differentiate  between  the  two  interface  conditions,  we  will  refer  to  (2.16)  as  a  pull  interface 
condition  and  (2.17)  as  a  push  interface  condition.  The  choice  between  the  two  equivalent  interface 
conditions  is  an  issue  of  implementation.  An  Eulerian  method,  such  as  the  finite-volume  method 
developed  in  Section  3.1.2,  combines  information  by  pulling  information  from  the  appropriate  bicar- 
acteristics  upwind  of  the  barrier.  A  Lagrangian  method,  such  as  the  particle  method  developed 
in  Section  4.2.3,  pushes  the  information  to  the  appropriate  bicharacteristics  located  downwind  of 
the  barrier.  Note  that  the  pull  interface  condition  (2.16)  is  a  many-to-one  function  and  the  push 
interface  condition  (2.17)  interface  condition  is  a  one-to-many  function.  In  contrast,  the  scattering 
for  the  classical  Liouville  equations  and  von  Neumann  equations  are  both  one-to-one  functions. 

We  assume  that  the  probability  of  a  particle  being  absorbed  by  the  barrier  is  zero  and  hence 
T(p)  +  R(p)  =  1.  By  defining 

1  if  \p(t~)\2  >  2m[V(x(t+))  —  V(x(t~))\  and 

T{p{t~))  = 

0  otherwise, 

i.e.,  total  transmission/reflection,  condition  (2.16)  reduces  to  the  classical  Liouville  condition  for 
which  bicharacteristics  are  uniquely  determined  for  each  ( x,p ).  When  T(p)  6  (0,1),  i.e.,  partial 
transmission/reflection,  the  bicharacteristics  are  no  longer  unique  and  instead  we  consider  multiple 
bicharacteristic  solutions. 

Every  interaction  with  a  barrier  potentially  introduces  a  reflected  and  transmitted  solution 
resulting  in  an  additional  bicharacteristic.  We  may  enumerate  the  solutions  and  define  a  bicharac¬ 
teristic  solution  to  the  Liouville  equation  as 

fk(x,p,t)  =  j  f(x,p)<pk(x,p,t;x,p)dxdp 

where 

<pk(x,p,t;x,p )  =  5(x(t)  -  x)6(p(t)  -p) 

is  the  fcth  global  bicharacteristic  for  H(x,p).  By  linearity  of  the  Liouville  equation  we  may  consider 
the  general  solution  as  the  superposition  of  the  bicharacteristic  solutions 

f(x,P,t)  =  ^ 2sk{H{x,p))fk(x,p,t ).  (2.18) 

k 

where  sk{H{x,p))  is  product  of  reflection  and  transmission  probabilities  along  the  Lth  bicharacter¬ 
istic. 
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Figure  2.1.  Particle  position  as  a  function  of  time  for  potential  V(x)  =  2\x\  —  H(x)  where  H(x)  is  the 
Heaviside  step  function.  Particle  has  initial  conditions  6(x  —  l)S(p). 


Except  for  simple  solutions  such  as  the  global-in-time  solution  for  a  piecewise-constant  poten¬ 
tial  or  the  local-in-time  solution  for  a  piecewise-quadratic  potential,  an  exact  solution  cannot  be 
explicitly  given.  Even  for  a  simple  discontinuous  oscillator  the  number  of  bicharacteristics  that 
need  to  be  tracked  becomes  cumbersome  in  a  short  time  interval.  See  Figure  2.2.  By  solving  the 
model  numerically,  we  mitigate  these  difficulties. 

2.3  Time  irreversibility  and  entropy 

The  Liouville  solution  at  time  t  for  an  initial  distribution  given  by  f(x,p,  0)  is  f(x,p,t )  = 

etC  f(x,p,  0)  where  the  Liouville  operator  C  =  V PH  •  Vx  —  V XH  ■  Vp.  By  reversing  the  momentum 

(p  i— >  —p),  we  have  f(x,—p,t )  =  e~tC  f(x,  —p,  0)  and  equivalently  f(x,—p,  0)  =  etcf(x,  —  p,  t). 

Hence,  the  Liouville  equation  is  time  reversible  under  the  negation  of  the  momentum.  Similarly, 

the  solution  of  the  Schrodinger  equation  at  time  t  for  an  initial  value  ip(x,  0)  is  ip(x,  t)  =  e~ltHip(x,  0) 

~  2 

where  the  Hamiltonian  operator  is  H  =  —  +  V(x).  By  reversing  the  phase,  i.e.,  by  taking  the 

complex  conjugate  of  the  probability  amplitude  (ip  i— >  ip),  we  have  ip(x,  0)  =  e~ltH ip(x,t).  Hence, 
the  Schrodinger  equation  is  time-reversible.  Unlike  the  Schrodinger  model  or  the  classical  Liouville 
model,  the  semiclassical  Liouville  model  does  not  preserve  time  reversibility. 

Consider  a  statistical  ensemble  of  particles — f(x,p,  0)  for  the  Liouville  equation  and  ip(x,  0) 
for  the  Schrodinger  equation.  Suppose  that  over  a  some  time  interval  [0,  t]  this  ensemble  collides 
with  a  barrier  resulting  in  the  scattered  solutions,  f(x,p,t)  and  ip(x,t )  respectively,  consisting 
of  transmitted  and  reflected  parts.  Now  consider  the  transmitted  and  reflected  solutions  with 
the  momenta  and  phases  reversed  as  initial  conditions,  namely  f(x,—p,t )  and  ip(x,t).  The  time 
evolution  traces  backward  along  the  bicharacteristics,  striking  the  barrier  and  resulting  in  scattered 
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solutions.  The  classical  Liouville  and  von  Neumann  solutions  correspond  to  the  original  ensemble 
f(x,—p,  0)  for  and  0);  the  semiclassical  solution  does  not.  The  classical  Liouville  equation  is 
time  reversible  because  there  is  only  one  bicharacteristic  for  each  initial  condition.  The  Schrodinger 
equation  is  time  reversible  because  the  solution  is  coherent,  containing  additional  phase  information 
which  results  in  constructive  and  destructive  interference  across  the  barrier.  The  semiclassical 
model  is  not  time  reversible  because  it  has  neither  of  these  properties — the  solution  is  completely 
decoherent  and  the  bicharacteristics  are  not  unique. 

Entropy  is  a  concept  closely  connected  to  time  irreversiblity.  The  statistical  entropy  S  of  a 
system  is  defined  as 

S{t)  =  £  Pi  log  Pi.  (2.19) 

i 

where  pi  is  the  probability  of  being  on  the  zth  local  bicharacteristc  at  time  t.  Whereas  the  clas¬ 
sical  Liouville  and  the  von  Neumann  equations  are  isentropic,  the  semiclassical  model  is  not.  As 
previously  stated,  the  scattering  relation  for  the  classical  Liouville  equations  and  von  Neumann 
equations  are  both  one-to-one  functions,  but  the  scattering  relation  for  the  semiclassical  model 
is  not.  The  semiclassical  interface  condition  mixes  information  from  different  bicharacteristics, 
thereby  increasing  the  entropy  of  the  system.  The  interface  condition  may  be  equally  posed  in 
terms  of  a  semiclassical  scattering  matrix 

T  R 
R  T 

which  relates  incident  and  scattered  states.  The  inverse  semiclassical  scattering  matrix 

T  R  \ 

T-R  R-T  I 

R  T 

R-T  T-R ) 

which  solves  the  time-reversed  semiclassical  scattering  problem  will  be  important  in  the  construction 
of  ghost  densities  for  the  slope  limiters  in  the  finite-volume  method  in  Section  3.1.2. 

As  an  example,  consider  the  harmonic  oscillator  with  a  delta-function  barrier 

V  ( x )  =  \x2  +  eaS(x) 

where  e  is  the  scaled  Planck  constant  and  a  is  some  parameter.  The  delta-function  quantum  barrier 
transmits  a  particle  with  momentum  p  with  probability  [1  +  (a/p)2]-1  [42].  Furthermore,  for  the 
Schrodinger  equation,  the  phase  of  a  transmitted  wavepacket  are  shifted  by  9t  =  —  tan-1  (a/p)  and 
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the  phase  of  a  reflected  wavepacket  is  shifted  by  Or  =  t  +  #t-  The  phase  is  also  shifted  by  it  when 
the  wavepacket  changes  direction  near  ±1.  Take  the  initial  conditions,  f(x,p,  0)  =  S(x  —  l)5(p). 
Then  the  particle  strikes  the  barrier  with  momentum  p  =  1  at  every  time  t  =  (n+  ^)vr  for  integer  n. 
Let  a  =  \/3.  Then  the  phase  shift  is  — 7t/3  for  transmission  and  ir/6  for  reflection.  A  semiclassical 
particle  with  momentum  p  =  1  will  be  transmitted  with  probability  j  and  reflected  with  probability 
Consider  the  probability  of  finding  the  particle  in  1R+  at  t  =  0,  7 r,  27t,  ....  For  the  semiclassical 
model,  it  is 

1  3  5  _9_  17  33  _65_ 

4’  8’  16’  32’  64’  128’ 

because  we  take  each  scattering  event  independently.  For  the  Schrodinger  equation,  it  is 

1  3  1  o  ±  3  x 

because  the  probability  amplitude  constructively  and  destructively  interferes  with  itself.  The  en¬ 
tropy  of  the  semiclassical  system  may  be  calculated  using  equation  (2.19).  We  have  that  at  t  =  nir, 
the  entropy  is 


0,  0.562,  0.661,  0.685,  0.691,  0.6926,  0.6930,  ... 

which  is  a  monotonically  increasing  function  asymptotic  to  the  equilibrium  value  log  2.  By  time 
t  =  2n  the  semiclassical  Liouville  solution  disagrees  with  the  Schrodinger  solution.  This  discrepancy 
is  due  to  limitations  of  the  semiclassical  model.  Clearly,  this  example  violates  assumptions  2 
and  4  on  page  8  that  require  each  interaction  with  the  quantum  barrier  to  be  independent. 

While  the  semiclasssical  solution  is  nonisentropic,  this  does  not  invalidate  the  model.  Since 
particles  at  the  classical  scale  interact  with  environment,  a  decoherent,  time- irreversible  solution  is 
a  physically  realistic  solution.  However,  a  coherent  solution  is  necessary  inside  a  quantum  region 
for  mesoscopic  and  periodic  potentials.  We  shall  examine  such  a  model  in  Chapter  5. 
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Chapter  3 

Semiclassical  Model  and  Numerical  Method  in  One  Dimension 

3.1  A  semiclassical  approach 

When  the  quantum  barrier  is  sufficiently  narrow,  the  barrier  may  be  modeled  using  the  time- 
independent  Schrodinger  equation.  We  may  then  derive  the  transmission/reflection  probabilities 
for  the  interface  condition  (2.16)  by  considering  the  current  density.  The  interface  condition  is  used 
to  connect  two  classical  domains  modeled  by  the  classical  Liouville  equation  (2.12). 

We  consider  an  algorithm  consisting  of  an  initialization  routine  and  a  Liouville  solver: 

1.  During  initialization,  we  determine  the  stationary  states  at  the  barrier  by  solving  the  time- 
independent  Schrodinger  equation.  The  solutions  may  be  found  by  considering  the  barrier 
as  an  open  quantum  system  [2]  outside  of  which  the  potential  is  constant.  Typically,  this 
may  be  done  by  using  a  quantum  transmitting  boundary  method  [28],  a  spectral  projection 
method  [32],  or  a  transfer  matrix  method  [1,  23,  15].  With  this  solution,  we  compute  the 
scattering  information,  namely  the  transmission  and  reflection  coefficients. 

2.  Following  initialization,  we  solve  the  Liouville  equation  using  a  finite  volume  method.  As 
done  in  [20]  the  interface  condition  (2.16)  is  built  into  the  numerical  flux  in  a  framework 
called  the  Hamilton  preserving  scheme.  This  yields  a  numerical  scheme  for  which  the  stabil¬ 
ity  condition — the  CFL  condition — is  hyperbolic,  namely  At  =  0(Ax,Ap)  with  l°°  and  l 1 
stability.  See  [20]. 

This  approach  aims  at  capturing  the  weak  limit  of  the  Schrodinger  and  von  Neumann  equations 
as  e  — >  0,  without  solving  the  Schrodinger  or  von  Neumann  equations  over  the  entire  domain, 
but  rather  just  at  the  quantum  barrier  and  only  in  the  initialization  step.  We  now  discuss  the 
initialization  routine  and  the  finite  volume  routine  in  detail. 
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3.1.1  Routine  initialization 

We  use  the  transfer  matrix  method  because  it  is  robust  over  a  wide  range  of  momenta.  On 
the  quantum  scale  we  decompose  a  one-dimensional  barrier  into  a  sequence  of  step  potentials  over 
which  we  solve  the  time-independent  Schrodinger  equation  exactly.  Take  a  quantum  barrier  in  the 
bounded  region  Q  =  [x\:x2]  and  take  the  potential  to  be  constant  outside  this  barrier — V{x)  =  V\ 
in  C\  =  (— oo,aq)  and  V{x)  =  V2  in  C\  =  (x2,oo). 

For  a  state  E  =  p2 /2m  the  time- independent  Schrodinger  equation 

—e2,ip"(x)  +  2  mV  (x)i/(x)  =  p2i/>{x) 


has  the  solution 


i/(x)  =  ^ 


aieiK i(x-xt)/e  +  &ig-iKi(x-x1)/e)  x  €  Cl> 

i/>Q,  X  G  Q 

a2e™ 2(x-x2)/e  +  &2g -iK2(x-X2)/e,  x  e  C2 


(3.1) 


where  =  \Jv2  ~  2mVi)2  and  the  coefficients  a\,  a2,  b\  and  b2  are  uniquely  determined  by 
the  boundary  conditions  at  X\  and  x2.  By  requiring  that  the  solution  i/{x)  and  its  derivative  be 
continuous,  ipQ  is  uniquely  determined  by  the  values  a\  and  b\  using  the  boundary  conditions 
iPq(x i)  and  iPq(x i).  In  turn,  the  values  a2  and  b2  are  uniquely  determined  by  the  values  i/q(x 2) 
and  iPq(x2).  Since  the  Schrodinger  equation  is  linear,  a2  and  b2  may  be  expressed  as  linear  functions 
of  a\  and  b\.  Hence,  for  each  momentum  p  we  may  relate  the  solution  in  C2  with  the  solution  C\ 
in  terms  of  the  transfer  matrix  M 


=  M 


m  n  m  12  /  a  1 

m2i  m2 2/  \bi 


(3.2) 


An  arbitrary  quantum  barrier  may  be  discretized  and  approximated  by  a  series  of  step  potentials, 
for  each  of  which  a  transfer  matrix  may  be  computed  analytically.  Specifically,  the  transfer  matrix 
may  be  approximated  as  M  =  Mn  •  •  •  M2Mi  with  Mj  =  D  '^2!  PjD'  '  2  where 

p  _  1  |  1  +  Kj/ 1  1  —  Kj/Kj+I 

\  1  —  Kj/Kj+1  1  +  Kj/Kj+1 
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Figure  3.1.  Approximation  of  a  potential  barrier  by  a  series  of  step  potentials.  The  effective  transfer 
matrix  M  =  Mn---M2Mi  where  Mj  is  the  transfer  matrix  for  a  step  potential  at  x:r 


is  the  transfer  matrix  associated  with  a  potential  jump  V(xf)  —  V(x-  )  and 


D;  = 


exp(iAxKj  /e) 


0 


0  exp(— iAxKj/s) 

is  the  transfer  matrix  associated  with  the  displacement  Ax  =  Xj  —  xj-\. 


(3.4) 


One  may  also  express  the  solutions  in  C\  and  C2  in  terms  of  a  scattering  matrix  S  which  relates 
the  incident  and  scattered  waves 


=  S 


-m2i/m22  l/ru-22 

A/m22  rni2/m22 


(3.5) 


J(x)  =  ^ 


(3.6) 


where  A  =  det  M  =  m22w.11  —  m^m^i-  By  considering  the  time  evolution  of  the  position  density 
p(x,t )  =  \il>(x,  t)\2  in  the  Schrodinger  equation,  one  derives  the  continuity  equation 

kp  +V'J  =  0 

where  the  current-density  is  defined  as  J(x)  =  em~l Im  (i/j'V'ip).  From  equations  (3.1),  one  has  that 

Ki  (|ai|2  —  |6i|2) /m,  x  €  Ci 

k2  (|«2|2  -  I&2I2)  /m,  x<eC2 

where  m  is  the  effective  particle  mass.  The  positive- valued  terms  of  the  J(x)  express  the  flux 
of  right-traveling  waves  and  the  negative-valued  terms  express  the  flux  of  left-traveling  waves.  In 
particular,  for  a  wave  incident  on  the  barrier  from  the  left  (62  =  0),  we  have  a2  =  tio-i  and  bi  =  r\a\. 
It  follows  that  the  reflection  coefficient  R\,  the  ratio  of  the  reflected  to  incident  current  densities, 
and  the  transmission  coefficient  Tj,  the  ratio  of  the  transmitted  to  incident  current  densities,  are 

R\  =  Ini2 


and  Tj  =  (a«2  /  ^1 )  |  ^1 1 2 


(3.7) 
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Similarly,  for  a  wave  incident  from  the  right 

R2  =  |r2|2  and  T2  =  (ki  /  K2)\t2\2 .  (3.8) 

The  transmission  and  reflection  coefficients  are  uniquely  determined  along  a  bicharacteristic.  It  is 
clear  by  time-reversibility  that  the  transmission  coefficient  along  any  bicharacteristic  is  independent 
of  direction 

Ti(p)  =  T2  (-y/p>  +  2m(V2-V1)')  .  (3.9) 

3.1.2  A  second-order  finite- volume  Liouville  solver 

Without  loss  of  generality,  we  shall  take  the  mass  m  =  1  in  which  case  we  equate  the  velocity 
with  the  momentum  p.  To  solve  the  semiclassical  Liouville  equation  (2.12),  we  use  a  Hamiltonian- 
preserving  finite- volume  method  [21].  We  consider  a  uniform  mesh  in  phase  space  with  grid  points 
at  (xi+1/2lpj+1/2)  and  denote  grid  spacing  Ax  =  xi+1/2  -  x{_ L/2  and  A p  =  pj+ 1/2  -  Pj-i/2  with 
i,j  E  Z.  Let  the  cell  centers  be  x.-L  =  \{xi+1/2  +  *j_1/2)  and  pj  =  \{Pj+ 1/2  +  Pj- 1/2).  For 
convenience  of  notation,  we  shall  take  po  =  0  and  p-j  =  —pj.  We  shall  consider  the  quantum 
barrier  to  be  located  at  a  cell  interface  xz+i/2  for  some  integer(s)  Z. 

Define  the  cell  average  over  the  cell  Cij  =  [®i_i/2, *i+i/2)  x  [pj-i/2,Pj+i/2)  as 

fS=A^JfCvf{x’p’tn)dxdp- 

The  hnite- volume  discretization  of  the  one-dimensional  Liouville  equation  (2.12)  is 

=  %  -  At  ppm'  -  dxVidpf%]  (3.io) 

where  the  discrete  operators  dxfij,  dpfij  and  dxVj  are 

dxfij  =  {f~+1/2,j  ~  ft i/2j)/Axi 
dpfij  =  (fi,j+ 1/2  -  /i,j-i/2)/Ap,  and 
Wi  =  (V-1/2  -  V+1/2)/Ax 
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with 


f± 

Ji+l/2,j 
fi,j  + 1/2 


rPj+ 1/2 


lim  5T-  /  f(x,p)dp, 


i+ 1/2 


Ap 


'Pj-1/2 


1  rxi+ 1/2 

—  /  f(x,pj+1/2)dx,  and 

JXi-1/2 


V 


i+ 1/2 


=  lim  V(a;). 


"i+1/2 


Upwinding  is  used  to  approximate  the  fluxes  /i^_1/2  7-  and  fi,j+ 1/2-  If  the  potential  V(x)  is  con¬ 
tinuous  at  some  point  a?j+i/2)  then  p(t+)  =  p(f_)  and  hence  f^+-l/2  3  =  /j+1/2  ^  which  reduces 
the  discretized  Liouville  equation  (3.10)  to  the  usual  upwind  finite  volume  scheme.  At  the  barrier 
xz+ 1/2  special  consideration  must  be  taken. 

From  conservation  of  the  Hamiltonian  (2.15)  we  have  that  the  incident  velocity  q3  (upwind  of 
the  barrier)  for  a  particle  transmitted  with  velocity  pj  is 


Qj 


Pj\J ,f  +  ^(I^Z+l/2  Vz+i/2) /'Pj\Pj 


Similarly,  the  transmitted  velocity  (downwind  of  the  barrier)  for  a  particle  incident  with  velocity 
Pj  is  —q~j.  The  incident  velocity  for  a  particle  reflected  with  velocity  pj  is  simply  —pj.  Note  that, 
whereas  —p~j  =  p3,  in  general  —q~j  7^  q3-  Further  note  that  by  time  reversibility  T(q_3)  =  T(pj) 
and  R{q~j )  =  R(pj)- 

The  left  and  right  limits  of  the  probability  distribution  /  in  the  cells  immediately  downwind  of 
the  quantum  barrier  are  determined  by  the  interface  condition  (2.16) 


fz+ 1/2,3  =  R^fz+l/2,-j  +  T(Qj)f(XZ+l/2’<lj)  for  3  >  0 

fz+l/2,j  =  R{Qj)fz+l/2,-j  +  T(Qj)f(xz+1/2^j)  for  3  <  0. 

The  values  for  f(x^+1j2,qj)  are  approximated  in  a  manner  similar  to  Scheme  II  of  [21].  Consider 
the  flux  incident  from  the  left  (q3  >  0) — the  same  treatment  applies  to  flux  incident  from  the  right. 
We  define  /(a^+1/2’ as  the  cell  average 

l  rij+i/2 

f(xz+i/2^j)  =  Pf(xz+i/2’P)dP  (3-n) 


Ph  1/2 


+  2(14 


z+1/2  ^z+1/2) 


where 
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The  integral  is  approximated  by  a  composite  mid-point  rule.  Since  the  limits  of  the  integral  are  not 
generally  gridpoints  in  the  p-direction,  some  care  must  be  taken.  If  Pk-1/2  <  Qj- 1/2  <  Qj+1/2  < 
Pk+1/2  for  some  k,  then  we  take 


f(xz+ 1/2’  lb')  /. Z+l/2,k  +  <ljcrp(f z+ 1/2, fc) 


where  the  slope  ctp(-)  in  the  p-direction  is  calculated  using  the  van  Leer  limiter 


ap(fij) 


f  -  fi,j- 1  \  ,  (  fij+l  -  ./' 


A  p 


(3.12) 


with  <j)(d)  =  (0  +  |0|)/(  1  +  |0|)  [29].  Otherwise  pfc_i/2  <  gj-1/2  <  ■  <  ^+1/2  <  Pfc+s+1/2  for  some 

k  and  s,  and  we  take 


f(xz+ 1/2’ 

1 


PjAp 


{(Pfc+1/2  -  Qj  — 1/2)  PkfZ+l/2,k  +  5  (Pfc+1/2  +  Qj-l/2)^p(Pkfz+1/2tk) 


+  Pk+l^Pf Z+i/2,k+l  ■*"■■■  “*"l,fc+s-1^^Z-|-l/2,fc+s-l 
+  (9j+1/2  -Pfe+s-l/2)  Pk+sf Z+i/2, k+s  +  i(Pfc+s-l/2  +  9j+l/2)<7'p(Pfe+s/z+l/2,fc+s)  }'  (3'13) 

For  a  second-order  accurate  method  we  use  a  slope-limited  piecewise-linear  interpolant  to  ap¬ 
proximate  the  right  and  left  density  limits 


/iipl/2 ,j  ~  fij  S1  2(1  Aj)A,TC7x(/ij) 


(3.14) 


where  A j  =  \vj\At/ Ax  and  the  slope  ax ( ■ )  in  the  ^-direction  is  calculated  using  the  van  Leer  limiter 


fij  -  fi- 


1,3 


Ax 


i  +  l,j 


fi 


fij  f i—l,j 


(3.15) 


Since  the  slope  <rx(-)  is  a  function  of  fi-ij,  fij  and  fi+ij  and  the  density  /  is  not  necessarily 
continuous  across  the  barrier  in  the  x-direction,  we  can  not  directly  use  (3.14)  and  (3.15)  to  calculate 
the  density  limits  at  the  barrier  interface.  Rather,  we  first  need  to  construct  the  ghost  densities  fz 
and  fz+i  across  the  barrier  using  the  scattered  densities  at  xz  and  x^+i  based  on  conservation  of 
mass.  Specifically,  downwind  of  the  barrier 


/z+1/2  =  /z+1/2(/z>  fz+iJz+2)  and  fz+1/2  =  fz+1/2(fz-i,fz,  fz+i) 
with  ghost  densities  fz  and  fz+1  located  upwind  of  the  barrier;  and  upwind  of  the  barrier 

fz- 1/2  =  ft-i/2(fz-ijz,  fz+ 1)  and  /“+3/2  =  fz+3/2(fz,  fz+iJz+2) 
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with  ghost  densities  /£  and  f%_ ^  located  downwind  of  the  barrier. 

Construction  of  the  ghost  densities  is  analogous  to  using  ghost  cells  to  enforce  inflow  and 
outflow  semipermeable  reflecting  boundary  conditions.  To  calculate  the  ghost  densities  upwind  of 
the  barrier  we  use  the  interface  condition  (2.16)  to  mix  together  the  densities  upwind  of  the  barrier 
that  will  subsequently  be  combined  through  transmission  and  reflection.  In  this  case 


fz,j  =  R(qj)fz+i,-j  +  T(qj)f(xz,  qj)  for  j  >  0, 

fz+ij  =  R(dj)fz-j  +  T(qj)f(xz+i,  qj)  for  j  <  0. 


To  calculate  the  ghost  densities  downwind  of  the  barrier  we  unmix  the  densities  downwind  of  the 
barrier  that  were  previously  combined  through  transmission  and  reflection  at  the  barrier.  In  this 
case 


T(Pj)f(xz+ 1,  ~ q-j )  ~  R(pj)fz,-j 
T(pj)  ~  Rip  j ) 

T{pj)f{xz,  -q-j)  -  R{Pj)fz+i,-j 
TiP 3 )  ~  RiPj) 


for  j  >  0, 


for  j  <  0. 


The  densities  f(xz+i,  ±q±j)  and  f(xz,Aq±j)  are  approximated  in  a  manner  similar  to  defini¬ 
tion  (3.11). 

To  approximate  f^+l/2  t°  second-order  in  the  p-direction  we  have 


± 

i,3z Fl/2 


fi,j  =F  li1  -  ^i)APapifij) 


with  A i  =  \dxVi\/S.t / /S.p  and  the  slope  crp(-)  defined  using  the  van  Leer  limiter  (3.12). 


3.2  Numerical  examples 

In  this  section  we  present  a  few  examples  of  both  pure  state  dynamics  and  mixed  state  dynamics 
in  order  to  verify  and  validate  the  semiclassical  model  and  numerical  scheme. 

For  a  mixed  state  solution  with  a  macroscopic  distribution,  we  are  not  limited  by  the  support  of 
the  wavepacket,  and  the  complexity  of  the  scheme  is  0((A.xApAt)_1)  where  Ax,  A p,  and  At  £■ 
For  direct  simulation  of  the  von  Neumann  equation,  not  only  must  we  resolve  e  in  space  and  time 
but  we  must  solve  the  equation  over  two  space  dimensions  and  one  time  dimension  so  the  complexity 
of  the  scheme  is  0(e-3).  When  e  <C  1,  the  computing  time  for  a  direct  von  Neumann  solution  is 
considerably  longer  than  for  the  multiscale  semiclassical  Liouville  solution. 
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The  numerical  Schrodinger  solution  may  be  computed  using  the  Crank-Nicolson  operator 

+  At)  =  (1  +  ie~x  AtH d)-1  (1  —  is-1  AtHr>)tl^(xi,t)  (3.16) 


where  the  discrete  Hamiltonian  operator 

— £2  1  —  2(5jj  +  Si ; 


TTd  = 


i+i 


+  V(Xi) 


(3.17) 


2m  (Ax)2 

with  Kronecker  delta  Su  =  1  and  Sij  =  0  if  i  ^  j.  Markowich,  Pietra  and  Pohl  [33]  show  that  for 
such  a  scheme,  in  order  to  guarantee  correct  approximation  to  observables  for  small  e.  one  needs  to 
take  Ax  =  o(s)  and  At  =  o(s).  One  may  also  compute  the  numerical  Schrodinger  solution  using  a 
pseudospectral  method  with  Strang  splitting  [4],  In  this  case,  one  splits  the  kinetic  and  potential 
terms,  so  that  for  each  time  step 

ip(x,t  +  At.)  =  eAtB/2J7~1 


eAtAT  eAtB/2^t) 


where 

A  =  -^—k2  and  B  =  —V(x) 

2  mi  is 

and  the  operators  T  and  T~x  denote  the  one-dimensional  discrete  Fourier  transform  and  discrete 
inverse  Fourier  transform  with  respect  to  the  x  and  k  variables.  One  can  use  a  mesh  that  is 
coarser  than  the  mesh  required  by  a  finite-difference  method  to  resolve  s  and  capture  the  correct 
dynamics  [4,  5].  Based  on  numerical  observation,  we  find  that  we  require  Ax  <  e/4  to  ensure 
numerical  convergence  to  the  correct  physical  observables  and  for  numerical  error  to  be  insignificant. 
When  the  potential  is  discontinuous,  we  find  that  the  solution  exhibits  artificial  oscillations  unless 
At  <  (A x)2/e  and  At  <  s/V(x). 

The  von  Neumann  equation 

ie-^p  =  Hp-  ( HpT)T  with  H  =  dxx  +  V (x) 
has  the  formal  solution 

p(x,  x'  ,t  +  At)  =  eieAt^  p(x,  x',  t)e~ieAtI{ . 


By  using  the  discrete  Hamiltonian  operator  (3.17),  we  may  approximate  the  von  Neumann  solution 
in  terms  of  the  Crank-Nicolson  operator  (3.16)  to  get  a  method  without  splitting  error 

/P+1  =  (1  +  is~1AtHD)~1(l  -  is~1AtHD)p*i  with 
p*j  =  (1  -  is~1AtHD)~1(l  +  is~1AtHD)p'lJi 
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where  p”-  =  p(xi,x'j,tn).  We  may  also  solve  the  von  Neumann  equation  using  a  pseudospectral 
method  with  Strang  splitting  [17], 


pn+l  =  eA  tB 


/2JF_1  ( eAtAT  (eAtB/2pn)) 


where 

A  =  -^—(k2  —  k'2)  and  B  =  —(V(x)  —  V(  x')) 

2  mi  ie 

and  the  operators  T  and  J-~x  denote  the  two-dimensional  discrete  Fourier  transform  and  discrete 
inverse  Fourier  transform  with  respect  to  the  ( x,x ')  and  ( k,k ')  variables.  The  FFTs  may  be 
optimized  by  exploiting  the  hermicity  of  the  density  matrix. 

Alternatively,  we  may  calculate  the  von  Neumann  solution  indirectly  by  solving  the  Schro- 
dinger  equation  for  several  states  and  then  using  definition  (2.7)  to  construct  the  density  matrix. 
This  simplifies  a  two-dimensional  problem  over  N2  gridpoints  to  n  independent  one-dimensional 
problems  over  N  gridpoints.  If  the  initial  distribution  is  localized  in  phase  space,  n  may  be  chosen 
to  be  appreciably  smaller  than  N,  saving  not  only  memory  but  also  contributing  to  a  considerable 
reduction  in  computation  time.  Furthermore,  this  approach  allows  us  to  implement  the  solution 
using  a  parallel  computer  cluster.  One  way  to  implement  such  a  scheme  is  to  use  states  generated 
by  taking  thin  slices  of  the  initial  distribution  along  the  ^-direction.  Consider  the  WKB  initial 
condition 

tp(x,0;x,p)  =  ((JxVZtt)-1/4  exp(— (x  —  x)2/Ag2)  ex-p(ipx/s),  (3.18) 

which  describes  a  wave  packet  with  an  0(1)  spread  in  position  and  0(e)  spread  in  momentum.  Let 
the  weight  distribution  in  the  definition  of  the  density  matrix  (2.7)  be 

f(x,p)  =  S(x  -  x0)  exp (— (p  -  p0)2/2SeO-2)/(s2(7pV/27r)  (3.19) 
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Using  the  Wigner  transform  (2.11),  we  have  the  equivalent  Liouville  initial  distribution 


f(x,p,0)  = 


1 


‘I'KG  x(Tp 


exp 


(x-Xq)2  ( P-PoY 


2al 


2 


(3.20) 


which  is  independent  of  e. 

To  compare  the  convergence  of  the  Schrodinger  and  von  Neumann  solution  in  the  semiclassical 
limit,  we  use  the  L1-error  of  the  position  probability  density  function  (pdf), 


\p(x,t)  —  \i/j(x,t)\2  dx 


with  p(x,t )  =  f^°  f(x,p,t)  dp.  We  replace  \ip(x,t)\2  with  p(x,x,t )  for  the  von  Neumann  solu¬ 
tion.  The  semiclassical  Liouville  model  should  also  predict  the  correct  weak  limit  for  multiphase 
solutions  when  interference  in  the  Schrodinger  and  von  Neumann  solutions  produce  oscillations  in 
the  probability  density  distribution.  To  measure  the  weak  convergence  in  the  semiclassical  limit, 
we  determine  the  .O-error  in  the  cumulative  distribution  function  (cdf),  i.e.,  the  antiderivative  of 
position  density  [14] 


n 


p(s,t )  -  \ip{s,t)\2  ds 


dx. 


In  each  example  we  compare  the  exact  or  numerical  semiclassical  Liouville  solution  with  nu¬ 
merical  Schrodinger  or  von  Neumann  solutions  for  equivalent  initial  distributions  and  potentials. 
Since  the  interactions  with  the  boundaries  are  not  relevant  to  the  study,  a  sufficiently  large  domain 
is  chosen  and  simulation  is  stopped  before  the  wave  envelope  reaches  the  boundaries. 


3.2.1  Schrodinger  0(1)  wave  envelope  with  a  step  potential 


Consider  the  step  potential 


V{x) 


0  if  x  <  0, 


(3.21) 


^  1  if  x  >  0. 

A  particle  impinging  on  this  potential  from  the  left  is  totally  reflected  when  the  incident  velocity 
is  less  than  1. 

We  find  the  exact  solution  by  the  method  of  characteristics  by  tracing  along  the  bicharacteristics 
backward  in  time  to  the  initial  conditions.  Let  Q(t)  =  {  (x,p)  \  x  <  0  and  x  —  pt  <  0,  or  x  > 
0  and  x  —  pt  >  0  }  be  the  region  in  phase  space  where  the  bicharacteristics  have  not  crossed  the 


quantum  barrier  at  x  =  0  within  a  time  t.  Then  the  exact  solution 
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f(x,P,t )  =  < 


f(x-pt,p,  0), 


(x,p)  G  fi(t) 


T  ■  f  I  -x  —  qt,  q,  0  +  R  ■  f  (—x  +  pt,  —p,  0) ,  otherwise 

p 


(3.22) 


where  the  incident  velocity  is  given  by  q  =  \Jp2  +  1  if  p  >  0  and  q  =  —  \/p2  —  1  if  p  <  0.  From 
equations  (3.3)  and  (3.7),  the  reflection  coefficient  is  given  by 


77  = 


P  +  9 


Note  that  when  p  G  [—1,0],  q  is  imaginary  and  R  =  1  indicating  total  reflection. 
Consider  the  WKB  initial  condition 


0)  =  ^4(x)ei5(x)/£ 


as  a  wavepacket  generalization  with  the  amplitude  and  phase  functions  given  by 

A(x)  =  (7 ra2)-1/4e-(— o)2/2a2 

S'(x)  =  ax2  +  bx  +  c. 

Since  S'(x)  is  a  quadratic  function,  we  can  calculate  the  Wigner  transform  of  ip(x,0)  exactly  to  get 

f{x,p,  0)  =  (IT £)-le(*-*o)2/*2e-(2ax+b-p)2/(e/a)2 

In  the  semiclassical  limit  (e  — >  0),  we  have 

f(x,p,  0)  =  A2(x)S{p  -  VxS(x)) 

=  (<7v^)_1e“(*_*o)a/<ra<5(p  -  (2ax  +  6)).  (3.23) 

By  taking  cr  =  0(1)  in  A(x),  we  create  a  wave  envelope  that  is  independent  of  s,  allowing  us  to 
study  the  convergence  of  solutions  as  e  — >  0.  When  a  ^  0,  the  distribution  of  phases  included  in 
the  Schrodinger  solution  is  also  0(1). 

Using  the  above  semiclassical  WKB  initial  conditions  (3.23),  we  note  that  when  t  =  —1/2 a 
the  position  density  for  the  Liouville  solution  (3.22)  exhibits  a  caustic  with  all  bicharacteristics 
intersecting  at  either  x  =  b/2a,  or  x  =  —b/2a  for  reflected  solutions.  Because  of  the  nonlinear 
change  to  the  incident  velocities,  the  transmitted  bicharacteristics  do  not  cross  simultaneously 
resulting  in  a  traveling  front,  the  leading  edge  of  which  is  unbounded. 
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Take  (xo,Po)  =  (—  1)  and  take  a  =  —  b  =  po  —  2a  =  |  and  cr  =  V.  Then  we  have  the  initial 

conditions 

^(®,0)  =  (io/7r)1/4e“100(x“xo)2ei(ax2+(po“2axo)x)£ 

for  the  Schrodinger  equation  and 

f(x,p,  0)  =  (10/7r)1/2e“200(x_xo)25(p-p0  -  2a(x  -  x0)) 

for  the  semiclassical  Liouville  equation.  The  numerical  Schrodinger  solution  is  solved  using  a 
Crank-Nicolson  finite-difference  method  over  the  domain  [— 1, 1]  with  mesh  size  Ax  =  At  =  10-6. 
The  exact  semiclassical  Liouville  solution  is  solved  by  tracking  characteristics  forward  in  time  with 
values  determined  by  the  initial  velocity  given  by  VS  =  |  —  \x.  We  compute  the  solution  at  time 
t  =  0.8. 

The  position  densities  for  several  values  of  e  are  shown  in  Figure  3.2.  The  convergence  results 
of  the  errors  in  the  two  solutions  are  listed  in  Table  3.1.  Based  on  this  study,  we  find  that  the  l 1 
convergence  rate  in  e  of  the  pdf  is  about  0.6  and  the  l 1  convergence  rate  in  e  of  the  cdf  is  about 
1.1. 

3.1.  Errors  in  solutions  of  Example  3.2.1  for  different  values  of  e. 


£ 

200-1 

800-1 

3200-1 

12800-1 

L-error  (pdf) 
d-error  (cdf) 

8.78  x  iCT1 
5.15  x  1(T2 

3.37  x  1CT1 
1.00  x  10-2 

1.55  x  10"1 
2.28  x  10-3 

8.61  x  10~2 
1.08  x  10~4 

3.2.2  Von  Neumann  solution  with  step  potential 

We  now  consider  the  solution  to  the  von  Neumann  equation  with  the  step  potential  given 
Example  3.2.1.  To  construct  a  von  Neumann  initial  condition  p(x,x',0)  which  corresponds  to  a 
Liouville  initial  condition  f(x,p,  0),  we  may  directly  use  the  definition  of  the  density  matrix  (2.7) 
for  some  weight  function  with  the  probability  amplitudes  ip(x,  t )  given  by  Gaussian  e-wavepackets 

ip(x,  0)  =  (7r£)-1/4e-(x-xo )2/^eiP0x/s  (3.24) 

The  Liouville  initial  condition  may  subsequently  be  calculated  by  a  Wigner  transform  of  the  density 
matrix.  Alternatively,  we  may  construct  the  density  matrix  by  using  the  inverse  Wigner  transform 
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applied  to  the  Liouville  initial  conditions  f(x,p,  0)  to  get 

/OO 

/(|(x  +  x'),p,  0)eip(x~x  M e  dp. 

-OO 

By  taking  the  Liouville  initial  conditions  to  be  the  Gaussian 


f(x,p,  0)  = 


1 


2tT  <7 x® p 


exp 


~(x  -  x0)5 

2cj2 


exp 


-(p-Pof 

2 


we  may  compute  the  von  Neumann  initial  conditions  exactly  to  get 

1  (  (\(x  +  x’)  -  xq)2  (x  -  x')2  ipo(x  -  x') 


p(x,x%  0)  = 


a 


,V2 7T 


exp 


2al 


2eVd2 


(3.25) 


(3.26) 


We  chose  ax  =  ap  =  0.05,  xo  =  —0.5  and  po  =  1.0  and  compared  the  solutions  to  the  von  Neu¬ 
mann  and  semiclassical  Liouville  equations  at  time  t  =  1.0.  The  von  Neumann  equation  was  solved 
using  the  psuedospectral  method  with  Strang  splitting  over  the  domain  [—1,1]  with  e  =  64“ 1, 
128-1,  256-1  and  512-1.  The  grid  spacing  was  fixed  at  Ax  =  2048-1  with  At  =  (Ax)2/e  to  ensure 
consistency  and  stability.  The  exact  solution  to  the  semiclassical  Liouville  model  is  calculated  using 
equation  (3.22). 

The  position  densities  for  the  semiclassical  Liouville  solution  and  the  von  Neumann  solution  for 
several  values  of  e  are  shown  in  Figure  3.3.  The  errors  in  the  two  solutions  are  listed  in  Table  3.2. 
Based  on  our  study,  we  find  the  convergence  rate  of  the  ^-error  of  the  pdf  is  about  0.7  as  e  — >  0 
and  the  convergence  of  the  ^-error  of  the  cdf  is  about  0.9  as  e  — >  0. 


3.2.  Errors  in  solutions  of  Example  3.2.2  for  different  values  of  e  . 


£ 

64-1 

I—1 

to 

OO 

1 

256-1 

512-1 

L-error  (pdf) 
d-error  (cdf) 

6.03  x  fO'1 
9.22  x  i0~2 

4.04  x  10"1 
4.83  x  10-2 

2.50  x  10”1 
2.53  x  10-2 

1.40  x  lO^1 
1.32  x  10~2 

We  may  also  consider  the  effect  of  incorporating  barrier  time  delay  in  the  approximation  of 
the  von  Neumann  equation  for  nonvanishing  e.  As  evident  from  the  offset  of  the  centers  of  the 
distributions  on  the  left  side  of  Figure  3.3,  one  source  of  error  is  the  time  delay  which  vanishes  in 
the  semiclassical  limit.  The  time  delay  may  be  considered  as  an  O(e)  correction  and  hence  we  may 
neglect  it  in  the  semiclassical  limit.  While  the  the  addition  of  a  delay  time  is  numerically  nontrivial, 
for  the  analytic  solution  (3.22)  it  is  a  straight-forward  modification. 
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Typically,  time  delay  is  defined  in  terms  of  the  Wigner  time  delay,  the  delay  to  the  group 

velocity  of  a  wave  packet  resulting  from  reflection  and  transmission.  As  such  it  is  meaningful  when 

the  wave  packet  has  a  well-defined  peak.  This  is  not  generally  the  case,  especially  when  the  barrier 

is  sufficiently  wide.  Considering  the  scattering  relation  (3.5),  the  reflection  and  transmission  group 

delay  times  for  unit  mass  are  [35] 

e  d  e  l  dt  e  d  e  1  dr 

rt  =  argi  =  -Im  )  and  rr  =  argr  =  -1m  ). 

p  dp  p  t  dp  p  dp  p  r  dp 

For  the  step  potential  (3.21),  we  have  from  equation  (3.3)  that  the  reflection  time  delay  is 

-ii 


Tr  =  2elm  [( pq )  ]  = 


pyjl  -  p 2 

when  p  E  [—1,0].  There  is  no  transmission  or  reflection  delay  time  for  p  ^  [—1,0].  To  incorporate 
the  time  delay,  we  make  the  replacement 


f(-x  +  pt,-p,  0)  -»•  f(-x+p(t  +  Tr),-p,  0) 

in  the  reflected  term  of  the  exact  solution  (3.22). 

We  compare  the  von  Neumann  solution  with  the  Liouville  solution  with  time  delay  correction. 
The  ^-errors  are  listed  in  Table  3.3.  Based  on  this  study,  we  find  that  the  addition  of  delay  time 
provides  some  improvement  to  the  model.  The  convergence  rate  of  the  / 1  -error  of  the  pdf  is  about 
1.3  and  convergence  rate  of  the  /Terror  of  the  cdf  is  about  0.9  as  e  — >  0. 


3.3.  Errors  in  solutions  of  Example  3.2.2  with  time  delay  correction. 


£ 

64-1 

I—1 

to 

00 

1 

256"1 

512-1 

d-error  (pdf) 
d-error  (cdf) 

3.67  x  KT1 
2.62  x  1(T2 

1.78  x  to-1 
1.65  x  1CT2 

7.05  x  10“2 
7.80  x  10“3 

2.23  x  10~2 
3.90  x  10~3 

3.2.3  Von  Neumann  solution  with  two  step  potentials 


We  may  consider  more  complicated  geometries  by  considering  multiple  barriers.  In  this  example 
we  construct  an  0(l)-wide  rectangular  barrier  by  taking  two  step  barriers  sequentially.  Consider 


V{x)  =  l 


\  if  x  E  [0,  |], 
0  otherwise. 
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We  take  the  initial  conditions  given  in  equations  (3.25)  and  (3.26)  with  ax  =  ap  =  0.05,  xq  =  —0.45 
and  po  =  1.1.  We  compute  over  the  domain  [—1.25, 1.25]  and  compare  the  solutions  at  time  t  =  1.2. 
The  von  Neumann  equation  is  solved  using  a  pseudospectral  method  with  Strang  splitting  as  in 
Example  3.2.2.  The  semiclassical  Liouville  solution  is  solved  using  the  numerical  method  proposed 
in  Section  3.1  using  N  grid  points  in  x  and  p  and  1.51V  steps  in  time.  The  results  are  shown 
in  Figure  3.4  with  e  =  0.002.  Even  with  a  fairly  coarse  mesh,  the  numerical  semiclassical  solution 
agrees  well  with  the  von  Neumann  equation  both  in  the  strong  limit  away  from  the  barrier  and  in 
the  weak  limit  between  the  two  step  potentials.  See  Figure  3.5. 

We  calculate  convergence  rate  as  Ax,  Ap,  At  — »  0  of  numerical  scheme  for  the  semiclassical 
Liouville  equation  by  computing  the  l1  -error  of  the  numerical  solutions  using  a  mesh  with  IV  =  50, 
100,  200,  and  400  grid  points.  For  an  “exact”  solution,  we  use  the  numerical  solution  using  N  = 
3200.  The  errors  are  listed  in  Table  3.4.  Based  on  this  study,  we  find  the  convergence  rate  of  the 
numerical  scheme  using  the  d-norrn  is  about  1.2. 

3.4.  Errors  in  solutions  of  Example  3.2.3  for  various  mesh  sizes  Ax. 


grid  points 

50 

100 

200 

400 

d-error 

3.32  x  10'1 

1.15  x  10'1 

4.72  x  10~2 

2.56  x  10-2 

3.2.4  Resonant  tunneling  von  Neumann  solution 

We  present  a  final  example  to  illustrate  a  specific  physical  model,  the  resonant  tunneling  diode 
(RTD)  [24,  34,  41],  An  RTD  consists  of  thin  layers  (a  few  nanometers  thick)  of  different  semi¬ 
conductors,  such  as  gallium  arsenide  (GaAs)  and  aluminum  gallium  arsenide  (AlGaAs),  that  are 
sandwiched  together  to  form  a  double-barrier  quantum  well  structure.  For  semiconductors  the 
de  Broglie  wavelength  is  on  the  order  of  tens  of  nanometers,  so  the  length  of  the  entire  RTD  struc¬ 
ture  is  on  the  length  scale  of  a  de  Broglie  wavelength.  The  region  outside  the  barrier  is  doped 
to  provide  a  sufficient  number  of  free  electrons.  Unlike  the  transmission  probabilities  of  the  step 
potentials  presented  the  previous  examples,  the  transmission  probability  of  an  RTD  is  not  a  mono¬ 
tonic  function  of  the  incident  particle  energy.  Rather,  it  is  oscillatory  and  admits  narrow  peaks 
of  total  or  almost  total  transmission  well  below  the  cutoff  energy  for  classical  transmission.  By 
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changing  the  bias  voltage  of  an  external  electrostatic  potential  applied  to  the  system,  the  resonance 
may  be  tuned  to  admit  electrons  of  varying  energies. 

We  shall  assume  that  the  electron  trajectory  is  ballistic.  In  the  quantum  region,  this  simpli¬ 
fication  is  appropriate  since  the  electron  mean  free  path  is  substantially  larger  than  the  barrier 
thickness.  However,  away  from  the  barrier  this  simplification  is  physically  unrealistic  since  the 
electron  mean  free  path  is  small  compared  to  the  classical  length  scale  for  a  dense  medium.  In  this 
case,  a  relaxation  term  or  collision  operator  should  be  added  to  the  Liouville  equation  to  capture 
the  particle  dynamics.  Since  we  require  that  the  Hamiltonian  be  only  locally  preserved,  the  model 
may  be  extended  to  a  dissipative  system,  for  which  the  Hamiltonian  is  continuous,  without  chang¬ 
ing  the  approach  discussed  in  Section  2.2  and  Section  3.1.  Hence,  for  the  purpose  of  validation,  the 
assumption  is  reasonable. 

We  construct  a  representative  barrier 

+5V0  x  £  (—00,  —a  —  b\ 

—  \Vqx/ (a  +  b)  +  14  x  £  (—a  —  b ,  —a]  U  (a,  a  +  b] 

—  |Vo x/(a  +  b)  x  £  [— a,  a  +  b] 

Vo  x  £  (a  +  6, 00) 

where  the  external  potential  bias  Vo  =  0.48,  the  thickness  of  each  barrier  b  =  0.9e,  the  thickness  of 
the  well  separating  the  barriers  2a  =  1.2e,  and  the  height  of  each  barrier  Vj,  =  2.25.  See  Figure  3.6. 
We  take  Gaussian  initial  distributions  (3.25)  and  (3.26)  with  ax  =  0.05,  ap  =  0.15,  xq  =  —  1  and 
Po  =  1.  The  solutions  are  computed  over  the  domain  [—4,4]  and  compared  at  time  t  =  2.5. 

The  von  Neumann  equation  is  solved  indirectly  using  the  WKB  initial  conditions  (3.18)  with 
weight  distribution  (3.19).  We  use  a  Crank-Nicolson  finite-difference  method  to  solve  the  Schro- 
dinger  equations.  To  ensure  that  the  weight  function  is  sufficiently  resolved,  we  take  N  =  (5e)_1 
Schrodinger  solutions  with  initial  values  equally  spaced  over  8ap  about  po- 

The  semiclassical  Liouville  solution  is  solved  using  the  numerical  method  proposed  in  Section  3.1 
using  an  N  grid  points  over  [—4, 4]  in  x,  2N  grid  points  over  [—3,  3]  in  p  and  3 N  steps  in  time. 
The  exact  solution  is  computed  using  equation  (3.22)  with  transmission  and  reflection  probabili¬ 
ties  calculated  using  the  transfer  matrix  method.  In  computing  the  transfer  matrix  for  both  the 
numerical  and  exact  solutions,  the  quantum  barrier  is  discretized  using  2000  grid  point  over  the 
length  6e  for  an  arbitrary  e.  The  results  are  shown  in  Figure  3.7  and  Table  3.5.  We  calculate  an 
^-convergence  rate  of  1.7  in  Ax,  Ap,  At. 
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3.5.  Errors  in  solutions  of  Example  3.2.4  for  various  mesh  sizes  Ax. 


grid  points 

80 

160 

320 

640 

T-error 

3.01  x  10_1 

1.37  x  10'1 

4.43  x  10~2 

8.90  x  10-3 

Figure  3.2.  Position  densities  for  the  semiclassical  Liouville  (top)  and  Schrodinger  (bottom)  solutions  of 
Example  3.2.1.  The  Schrodinger  solution  shows  e  =  (a)  200”1,  (b)  800-1,  (c)  3200-1  and  (d)  12800-1.  The 
position  density  of  Liouville  solution  exhibits  a  caustic  near  x  =  0.08  and  the  peak  is  unbounded.  For  the 
Schrodinger  solution  the  peak  reaches  a  height  of  19  for  the  e  =  12800-1. 
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Figure  3.3.  Position  densities  for  the  semiclassical  Liouville  (top)  and  von  Neumann  (bottom)  solutions  of 
Example  3.2.2.  The  von  Neumann  plot  shows  e  equal  to  (a)  64_1,  (b)  128-1,  (c)  256-1  and  (d)  512-1. 
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Figure  3.4.  Position  densities  for  the  numerical  semiclassical  Liouville  (top)  and  von  Neumann  (bottom) 
solutions  of  Example  3.2.3.  The  •  in  the  Liouville  plot  shows  the  numerical  solution  for  with  150  grid 
points  over  the  domain  [—1.25, 1.25].  The  solid  line  shows  the  numerical  solution  for  3200  grid  points.  The 
von  Neumann  solution  is  for  £  =  0.002. 


Figure  3.5.  Detail  of  Figure  3.4  showing  position  densities  for  the  numerical  semiclassical  Liouville  and 
von  Neumann  solutions.  The  •  shows  the  numerical  solution  for  with  150  grid  points  over  the  domain 
[—1.25,1.25].  The  solid  line  shows  the  “exact”  Liouville  solution  and  the  von  Neumann  solution  using 
£  =  0.002. 
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Figure  3.6.  Transmission  probability  as  a  function  of  the  momentum  p  for  the  RTD  barrier — shown  in  the 
inset — presented  in  Example  3.2.4. 


-3-2-101234 


x 


x 


Figure  3.7.  Position  densities  for  the  numerical  semiclassical  Liouville  (top)  and  von  Neumann  (bottom) 
solutions  of  Example  3.2.4.  The  Liouville  solution  shows  the  numerical  solution  for  (a)  80,  (b)  160, (c)  320 
and  (d)  640  grid  points.  The  von  Neumann  solution  is  plotted  for  e  =  (a)  50_1,  (b)  100-1,  (c)  200~1  and 
(d)  400-1.  The  exact  solution  (*)  is  also  plotted  in  each  case. 
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Chapter  4 

Semiclassical  Model  and  Numerical  Method  in  Two  Dimensions 

4.1  Extension  to  multiple  dimensions 

Multiple  dimensions  presents  several  challenges  to  computing  both  the  quantum  von  Neumann 
equation  and  the  semiclassical  limit.  Not  only  do  the  algorithms  become  more  complicated  as  more 
variables  must  be  taken  into  account  but  also  limitations  in  computer  resources  (computing  time 
and  memory)  introduce  unavoidable  obstacles.  Such  obstacles  are  the  primary  motivation  for  the 
development  of  a  computationally  efficient  and  tractable  semiclassical  model.  The  von  Neumann 
model  for  d-dimensional  dynamics  requires  a  2d-dimensional  density  matrix.  Whereas  one  may 
need  15MB  of  computer  memory  to  compute  one-dimensional  dynamics  over  a  unit  interval  with 
the  £  =  500-1  using  a  direct  method,  one  would  need  15TB  of  computer  memory  to  compute 
the  equivalent  two-dimensional  dynamics  using  the  same  method.  Aside  from  memory,  a  two- 
dimensional  solution  needs  a  million  times  as  many  floating  point  operations  as  the  equivalent 
one-dimensional  solution.  While  an  indirect  method  of  solution  mitigates  the  memory  concern  by 
solving  a  large  number  of  d-dimensional  Schrodinger  equations  independently,  such  an  approach  is 
impractical  for  general  initial  distributions. 

The  so-called  “curse  of  dimensionality”  also  afflicts  the  numerical  solution  to  the  semiclassical 
model.  One  could  solve  most  one-dimensional  problems  with  a  typical  computer  using  a  dense, 
concurrent  finite  volume  approach.  In  higher  dimensions,  such  an  approach  in  general  is  simply 
not  effective.  Consider  the  solution  to  a  two-dimensional  problem,  which  requires  four  dimensions 
in  phase  space.  Even  a  rather  coarse  mesh  using  an  array  of  1004  floating-point  numbers  requires 
380MB  of  memory.  Since  we  also  require  an  additional  “swap”  array,  we  find  that  1004  is  a  practical 
limit  for  brute  calculation.  Because  a  typical  problem  requires  at  least  100  grid-points  over  a  unit 
interval  to  resolve  details  (see  Figure  3.7),  the  finite-volume  method  developed  in  Section  3.1.2  is 
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ineffective  for  general  multi-dimensional  semiclassical  problems.  A  sparse  matrix  algorithm  may  al¬ 
leviate  some  of  the  difficulty;  however,  such  an  approach  is  viable  only  when  the  density  information 
is  sufficiently  local  (such  as  a  front),  which  is  typically  an  exception  for  von  Neumann  solutions.  In 
addition,  sparse  matrices  introduce  numerical  record-keeping  issues  further  reducing  the  numerical 
efficiency  of  the  approach. 

Instead,  we  consider  a  mesh-free,  particle  method  as  an  effective  alternative.  For  non-interacting 
particles,  the  bicharacteristic  solutions  may  be  computed  independently  thereby  eliminating  the 
memory  constraints.  While  a  finite-volume  approach  requires  a  concurrent  solution  using  a  dense 
array,  a  particle  method  algorithm  may  be  easily  adapted  for  parallel  computation  on  a  distributed 
computer  cluster  reducing  the  simulation  run  time.  Furthermore,  because  other  related  physical 
models,  such  as  for  plasmas,  often  rely  on  particle  methods  for  simulation  it  is  quite  natural  to  use 
such  an  approach  for  thin  quantum  barriers. 

While  mitigating  one  set  of  challenges,  the  particle  method  introduces  another  set.  Since  we 
use  the  bicharacteristics  to  track  information  directly,  divergence  of  the  bicharacteristics  is  prob¬ 
lematic  for  all  but  trivial  examples.  Because  of  this  one  must  periodically  reconstruct  of  the  data. 
Furthermore,  reconstruction  of  the  data  is  difficult  in  regions  where  the  particles  are  sparse  and 
smoothing  techniques  may  required  to  mollify  the  numerical  solution. 

Formally,  we  define  a  particle  as  the  approximation  to  a  Dirac  measure  using  some  type  of 
cutoff  function  [37].  The  particle  method  consists  of  first  approximating  the  initial  conditions 
/o(r)  =  f  fo(r)8(r  —  r)  dr  by  a  linear  combination  of  Dirac  measures  /g  =  1  Wj8{r  —  rj)  for 

some  set  {rj,Wj}  with  position  rj  £  Rd  x  M.d  and  weight  Wj  >  0  where  N  is  the  sample  size. 
The  set  {r3 ,  wj }  may  be  chosen  by  either  a  Monte  Carlo  method  or  a  deterministic  method.  In  a 
Monte  Carlo  method,  one  samples  r3  randomly  from  a  distribution  and  sets  Wj  =  N~x  f  fo(r)dr. 
In  a  deterministic  approach,  one  assigns  rj  based  on  a  uniform  or  nonuniform  mesh  and  sets  Wj  = 
fc  fo(r)  dr  for  a  cell  Cj  £  M.d  x  Md.  A  problem  is  solved  by  considering  the  time  evolution  of  these 
particles  with  the  appropriate  weights.  To  solve  the  Liouville  equation,  where  8(x(t.)  —  x)8(p(t)  —  p) 
defines  a  single  bicharacteristic  for  the  Hamiltonian  H(x,p),  we  solve  the  Hamiltonian  system  of 
equations  (2.1)  for  each  particle  sampled  from  f0(x,p). 

The  focus  of  this  chapter  is  to  develop  an  efficient  numerical  discretization  of  the  semiclassi¬ 
cal  model  presented  in  Section  2.2  in  two  dimensions.  Although  we  limit  the  discussion  to  two- 
dimensional  physical  space,  the  extension  three  dimensions  follows  using  a  similar  treatment.  As 
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in  the  one-dimensional  case,  we  use  the  simplifying  assumptions  that  the  effective  width  of  a  thin 
quantum  barrier  is  0(s),  the  distance  between  neighboring  barriers  is  0(1),  and  the  gradient  in 
the  the  potential  is  0(1)  except  at  quantum  barriers.  Furthermore,  we  require  that  coherence  time 
be  sufficiently  short  so  that  we  may  neglect  interference  away  from  the  barrier.  In  two-dimensional 
space,  we  consider  the  quantum  barrier  as  a  smooth  one-dimensional  curve  F q  separating  two  clas¬ 
sical  regions  C\  and  C?,.  In  addition  to  changing  across  the  width  of  the  barrier,  the  potential  may 
also  change  along  the  length  of  the  barrier  at  either  the  classical  0(1)  length  scale  or  a  quantum 
O(e)  length  scale.  Hence,  in  the  semiclassical  limit  not  only  is  the  potential  discontinuous  at  the 
barrier  in  a  direction  normal  to  the  barrier  curve,  but  the  potential  may  also  be  discontinuous  along 
the  barrier  curve.  As  in  the  one-dimensional  case,  we  prescribe  an  interface  condition  to  match  local 
solutions  in  order  to  construct  the  global  bicharacteristic  solution.  But  unlike  the  interface  condi¬ 
tion  for  one-dimensional  dynamics  that  combines  just  two  bicharacteristics,  the  interface  condition 
for  two-dimensional  dynamics,  potentially  joins  a  continuum  of  bicharacterist  computed  using  the 
time-independent  Schrodinger  equation. 


4.2  A  semiclassical  approach  and  numerical  discretization 
4.2.1  A  semiclassical  approach 


Whereas  in  one  dimension  for  which  there  are  only  two  momenta  associated  with  a  Hamil¬ 
tonian  E  (namely  +p  and  —p),  in  two  dimensions  there  are  a  continuum  of  momenta  p(0)  = 
(pcos91psin9).  Along  a  local  bicharacteristic  the  momentum  of  a  particle  is  uniquely  determined 
by  continuity  of  the  potential.  But  across  a  quantum  barrier,  where  potential  is  discontinuous 
and  the  gradient  of  the  potential  is  classically  undefined,  the  continuation  of  the  momenta  is  not 
unique.  In  order  to  match  bicharacteristics,  we  use  information  at  the  quantum  scale  to  construct 
an  interface  condition. 

The  two-dimensional  analogues  to  the  pull  interface  condition  (2.16)  and  the  push  interface 
condition  (2.17)  are 

/7t/2  /» 7t/2 

R(9-pin,6out)f{xin,pin,9)d9  +  /  T(0;  q-m,  0out)/(xin,  gin,  9)  d9  (4.1) 

-tt/2  J— tt/2 

and 


/7t/2  r7r/‘2 

#(0; Pout,  0in)/(xout, Pout,  0)  d6+  /  T(9-,qout,9in)f(xout,qOut,0)d6  (4.2) 

-7r/2  J  —  7t/2 


tt/2 

-tt/2 
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respectively.  Here,  R{0Out',  Pin,  #in)  is  the  probability  of  a  particle  incident  with  momentum  p-m  at 
an  incident  angle  0-iri  being  reflected  at  a  reflection  angle  0Out,  T(0out-,pm,0m)  is  the  probability  of 
a  particle  incident  with  momentum  p-m  at  incident  angle  0;n  being  transmitted  at  a  refraction  angle 
0out,  and  q2  =  p2  —  2 mAV. 

Where  the  potential  is  discontinuous,  one  may  treat  the  gradient  of  the  potential  as  an  impulse 
force.  However,  the  direction  of  such  a  force  may  not  be  well-defined  at  the  classical  scale.  If  the 
potential  is  discontinuous  both  in  the  direction  normal  to  the  barrier  and  also  along  the  length  of 
the  barrier,  we  must  use  the  solution  at  the  quantum  scale  to  determine  the  appropriate  scattering 
angles.  We  shall  refine  this  idea  when  we  discuss  the  quantum  scale  solution  in  Section  4.2.2.  If 
the  semiclassical  potential  V (x,  y )  is  discontinuous  in  the  direction  normal  to  the  quantum  barrier 
curve  rQ  but  is  continuous  along  the  length  of  Tq,  we  take  the  impulse  force  normal  to  the  barrier 
curve.  In  this  case,  one  has  as  a  consequence  of  conservation  of  the  Hamiltonian  that  the  change 
in  momentum  for  a  reflected  particle  is 


Ap  =  —  2(pin  •  n)n  (4.3) 

where  n  is  the  unit  normal  to  the  barrier  pjn  denotes  the  incident  momentum;  and  for  a  transmitted 
particle,  the  change  in  momentum  is 

Ap  =  ( \/ 1  Pin  •  n\2  +  2mAV  -  pin  •  n)n.  (4.4) 

One  may  relate  the  angle  of  refraction  to  the  angle  of  incidence  (defined  with  respect  to  the  unit 
normal)  by  using  the  conservation  of  the  Hamiltonian  to  derive  an  expression  analogous  to  Snell’s 
law  of  geometric  optics 

sin 02  =  sin 01  / y/l  +  2m(V1  -  V2)/\p>\2 

where  0\  is  the  angle  of  incidence,  02  is  the  angle  of  refraction,  V\  is  the  potential  on  the  incident 
side  and  V2  is  the  potential  on  the  scattered  side.  From  this  expression,  one  may  note  that  when 
the  angle  of  incidence  is  greater  (shallower)  than  a  critical  angle 

0i  >  Oc  =  cos-1  (V 2m(V2  -  Hi)/|p|2)  (4.5) 

the  particle  is  totally  reflected  by  the  barrier. 

In  the  following  sections,  we  present  the  particle  method  which  solves  the  semiclassical  Liouville 
equation.  As  in  the  one-dimensional  case,  the  algorithm  consists  of  an  initialization  routine  and 
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a  Liouville  solver.  During  initialization,  we  determine  transmission  and  reflection  coefficients  as  a 
function  of  the  incident  momentum  along  the  interface  from  both  sides.  Calculation  of  transmission 
and  reflection  coefficients  for  a  two-dimensional  interface  is  an  extension  of  calculation  of  those 
coefficients  in  the  one-dimensional  case.  As  in  one  dimension,  we  solve  the  problem  by  computing 
the  stationary  solution  to  the  time-independent  Schrodinger  equation.  For  the  semiclassical  model, 
we  consider  the  quantum  barrier  as  a  curve  Tq  separating  two  classical  regions  C\  and  C2.  Because 
the  potential  may  change  along  the  length  of  the  curve,  we  compute  the  transmission  and  reflection 
coefficients  locally  at  each  point  along  the  curve.  Consider  a  point  (xo,yo)  6  Tq  and  define  the 
local  coordinates  (, x,y )  where  the  ^-direction  is  normal  to  Tq  the  barrier  at  (xq ,  yo )  and  the  y- 
direction  is  parallel  to  Tg  at  (xo,yo).  By  assumption,  the  width  of  the  quantum  barrier  is  0(e)  in 
the  ^-direction  and  the  length  of  the  quantum  barrier  is  0(1)  in  the  y-direction.  Formally,  we  will 
associate  the  semiclassical  quantum  barrier  Tq  with  a  region  Q  bordered  by  the  classical  regions 
C\  and  C2.  By  assumption  the  gradient  of  the  potential  V(x,y)  in  classical  regions  is  0(1). 

Consider  the  two-dimensional  time-independent  Schrodinger  equation 

(Jj^2  +  ^£(x,y)  +Ve(x,y)'ip£(x,y)  =  E^£ (x ,  y) .  (4.6) 

By  rescaling  x  and  y  by  s  (x  =  ex  and  y  =  ey),  the  Schrodinger  equation  (4.6)  may  locally  be 
expressed  as 

1  /  <92  d2  \ 

”2^  V5i2  +  Q~2  J  ^y)  +  V(x,y)iJj(x,y)  =  E^(x,y).  (4.7) 

By  letting  e  —>  0,  we  may  regard  C\  and  C2  as  the  semiinfinite  regions  C\  =  {  (x,  y)  \  x  <  x\  } 
and  C2  =  {  (x,y)  \  x  >  X2  }  separated  by  an  infinite  strip  Q  =  {  (x,  y)  \  x\  <  x  <  X2  }  for 
some  x\  and  X2-  We  solve  the  time-independent  Schrodinger  equation  over  Q  using  information 
in  C\  and  C2  to  generate  transmission  and  reflection  coefficients.  We  are  interested  in  computing 
the  transmission  and  reflection  coefficients  locally,  so  variations  in  the  potential  that  are  on  the 
classical  0(1)  lengthscale  in  the  y-direction  may  be  neglected  at  the  quantum  O(e)  lengthscale. 
Hence,  we  define 


V(x,y) 


Vi,  (x,y)  €  Ci 

<  VQ(x,y),  (x,y)  €  Q 
V2,  ( x,y)eC2 


(4.8) 


where  Vi  and  V2  are  constants. 
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4.2.2  Routine  initialization 

In  this  section  we  discuss  the  quantum  transmitting  boundary  (QTB)  method  [7,  28]  as  a  means 
of  determining  the  reflection  and  transmission  coefficients  of  the  thin  two-dimensional  quantum  bar¬ 
rier.  The  QTB  method  is  used  to  solve  the  time-independent  Schrodinger  equation  in  a  region  with 
open  boundary  conditions.  By  using  continuity  of  the  solution  and  its  derivative  at  the  boundaries 
of  an  open  quantum  system  in  conjunction  with  a  solution  with  undetermined  coefficients  in  the 
exterior  region,  one  formulates  a  boundary  value  problem  for  the  interior  region.  The  unknown  co¬ 
efficients  are  eliminated  from  the  problem  by  combining  the  Dirichlet  boundary  conditions  with  the 
Neumann  boundary  conditions  to  get  mixed  boundary  conditions.  Once  the  solution  in  the  interior 
is  known,  it  may  be  used  on  boundaries  to  recover  the  unknown  coefficients.  The  approach  differs 
from  the  transfer  matrix  method  in  that  the  QTB  method  first  solves  the  Schrodinger  equation 
numerically  in  the  quantum  region  as  a  means  of  deriving  the  scattering  coefficients.  The  transfer 
matrix  method,  on  the  other  hand,  uses  the  exact  solution  to  the  Schrodinger  equation  over  several 
smaller  regions  as  a  means  of  connecting  boundary  values.  While  the  transfer  matrix  method  is 
accurate  over  a  wider  range  of  momenta,  the  QTB  has  a  natural  extension  to  two  or  more  dimen¬ 
sions.  Approaches  to  compute  the  solution  to  the  time-independent  Schrodinger  equation  using 
transfer-matrix  methods  for  two-dimensional  geometries  are  discussed  in  [10,  39]. 

We  adapt  the  approach  proposed  by  Lent  and  Kirkner  [28].  Consider  the  solution  to  the  local 
time-independent  Schrodinger  equation  (4.7).  Here,  and  in  the  sequel,  the  tildes  on  x  and  y  are 
dropped  in  order  to  simplify  notation.  Without  loss  of  generality,  we  take  the  potential  in  region 
Ci  to  be  zero  (Hi  =  0).  In  this  case,  E  =  p\/2m  =  p\/2m  +  Hz  where  p\  is  the  magnitude  of  the 
momentum  of  a  particle  in  region  Ci  and  P2  is  the  magnitude  of  the  momentum  of  a  particle  in 
region  C2. 

The  solution  to  the  local  Schrodinger  equation  (4.7)  may  be  written  as  the  piecewise  function 

Vh  (x,y),  (x,y)€Ci 
^(x,y)  =  lipQ(x,y),  ( x,y)eQ 
^2  (x,y),  (x,  y)  €  C2 

for  which  the  components  ipi (x,y),  ipQ(x,y)  and  ip2(x,y)  are  related  by  appropriate  matching 
conditions.  In  regions  Ci  and  C2,  where  the  potential  V(x,y)  is  constant,  the  Schrodinger  equation 
simplifies  to  the  Helmholtz  equations 


Aipj(x,  y)  =p2jipj(x,y), 


3  =  1,2 


(4.9) 
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which  have  the  general  solutions 

ipj{x,y)  = 


COB  9+yainO)  dQ ,  j  =  12 


(4.10) 


The  current  density  as  discussed  in  Section  3.1.1  is  defined  as  J(x,y)  =  m  1Im  ('il>(x,y)'V'i/j(x,y)') 
where  m  is  the  effective  mass.  For 


ip(x,  y)  =  f  a(6)eip{xcosd+ysin9)  d6 
J  — 7T 


we  have 


/77  /»7T 

a(0)e-ip(xcose+ysind)  d6  /  (cos  9,  sin  9)a(9)eip^x cos  e+y  sin e)  dd 

-77  J  —  77 

a(0i)a(02)(cos  02 ,  sin  02)e-ip(x(cos  dl  ~cos  e^+y(sin  6,1  _sin  dQldQ2_ 


’  —77 
f*77  (‘77 


—  ip 


'  —77  J  —77 


To  determine  transmission  and  reflection  coefficients,  we  need  the  average  directional  flux  at  the 
barrier.  Integrating  this  over  y  we  have 

/OO 

^{x,y)Vi){x,y)  dy 

-OO 

a(91)a{62)(cos02,sm92)e-ipx(cosei-cosd2)e-ipy{sin9l-sine2UydO1dh 
a(6»i)a(6»2)(cos02,sin6»2)e-ipx(cosei-cose2)(5(0i  -92)d91d92 
=  ip  I  |a(#)|2(cos  9,  sin#)  d6 

J  —77 

Hence,  the  directional  contribution  to  average  current  density  along  the  y-axis  is 


=  ip 


=  ip 


p77  (‘77  pOO 


'  —77  J  —77  J  —OO 
(‘77  p77 


'  —77  J  —77 


Jj(x,y,6)  =  \aj(9)\2pj (cos  6,  sin  9)  j  =  1,2 


(4.11) 


which  says  that  the  magnitude  of  the  particle  flux  through  a  point  in  region  C\  at  an  angle  9  is 
Pi  | a\  (0)| 2  and  the  magnitude  of  the  particle  flux  through  a  point  in  region  C2  at  an  angle  9  is 
p2\a2(9)\2. 

Consider  particle  initially  in  region  C\  that  strikes  the  barrier  from  the  left  with  momentum  p\ 
at  an  angle  of  incidence  #;n.  The  particle  scatters  with  momentum  pi  into  region  C\  if  reflected 
and  momentum  p2  into  region  C2  if  transmitted.  In  this  case,  the  solutions  to  equations  (4.9)  are 


^i(x,y) 

^2 (x,y) 


_  qWi({x—xi)  cos  #in+y  sin  $in) 


f7r/2 


'  —77/2 


r(9)e~ipi^x~xi^  cosd+y sin61)  cW 


f-K/2 


'  —77/2 


f  ^eip2((x-x2)  COS  e+y  sinB) 


(4.12a) 

(4.12b) 
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where  r(9)  and  t{9)  are  some  yet  unknown  scattering  distributions.  The  probability  that  a  particle 
is  scattered  at  some  angle  equals  the  ratio  of  the  scattered  current  density  to  the  incident  current 
density.  From  (4.12), 

/7r/2  r7r/‘2 

\r (9)\2 pi(cos  9,  sin  6)  d6  +  /  \t(9)\2p2(cos9,sin9)  d9. 

-7r/2  J—ir/2 

Hence,  the  reflection  and  transmission  probability  distributions  over  the  sector  ( 6  —  d d9,9  +  \d9) 
for  incident  (p\ ,  9\n)  are 

dR{9)  =  \r{6)\2^^-d6  and  dT{9)  =  \t{9)\2  p2COii®  d9. 

cos  9m  pi  cos  9-in 

While  the  form  of  the  solutions  (4.12)  is  convenient  for  discussing  scattering  solutions,  it  is 
inconvenient  for  actually  determining  them  since  the  unknows  r{9)  and  t{9)  are  coupled  through 
the  integrals.  The  Schrodinger  solution  in  region  Q  is  a  boundary  value  problem  with  boundaries 
parallel  to  the  y-axis.  By  expressing  the  solutions  (4.12)  in  terms  of  the  y-components  of  the 
momenta,  we  may  rewrite  them  in  the  equivalent  forms 


si{0e~immx~xi)e^v  d^, 

-OO 

/OO 

S2(£,ymmx-x*)^v  dt, 

-OO 


with  i  =  p\  sin  0 
with  £  =  p2  sin  6 


where  the  x-components  of  the  momenta  are 

m(0  =  \JpI-£2  and  V2(0  =  \jv\-i2-, 


(4.13a) 

(4.13b) 


the  complex  scattering  coefficients  are 


si(0 


r(9)picos9,  if  |£|  <  pi 
<  and  s2(£) 

0,  otherwise 


t(9)p2  cos  9,  if  |£|  <  p2 

< 

0,  otherwise; 


(4.14) 


and  y-component  of  momentum  of  the  incident  particle  is  £in  =  pisin0in.  Note  that  r{9)  = 
si(0?7i(0  and  =  s2  (0V2(0  for  77 1(£)  and  772  (6)  real.  By  taking  the  Fourier  transform  of  the 
solutions  (4.13)  with  respect  to  y  we  have 


teO  =  -  Ueim^){x~xi)  +  si(Oe 

Mx,C)  =  s2(£)eir72(0(x'X2)- 


-iv  dOix-x-t) 


(4.15a) 

(4.15b) 
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The  Fourier  transform  of  the  Schrodinger  equation  (4.7)  is 


—Z-2$q(xiQ  +  2m  /  VQ(x,y)ip(x,y)e  <y  dy  =  0. 


(4.16) 


where  we  define 


/OO 

^(x,y) 

-OO 


e~^y  dy 


to  be  the  Fourier  transform  of  ij)Q{x,y)  analogous  to  our  definitions  for  Vh(cc,£)  and  i>2(x,  £)• 

By  requiring  that  the  solution  ip(x,y )  and  its  first  derivatives  be  continuous,  we  have  the 
matching  conditions  at  x  =  x\  and  x  =  X2 


*j  A.  Q  *  d  ~ 

=  tpQ(xj,0  and  —ipj(xj,€)  = 

for  j  =  1,2.  Applying  these  matching  conditions  to  equations  (4.15)  we  have 


(4.17a) 


i>Q{x  1,6  =  ^-^in)  +  Sl(0) 


^q(x2,0  =  S2(0, 


r\  r\ 

-q^q{xi,0  =  irh(0s(t  -  £in)  -  *»7i(0si(0)  =  'im^)s2(0- 

Eliminating  the  unknowns  «i(£)  and  $2(£)  gives  the  boundary  conditions 


(4.18a) 

(4.18b) 


imiOi’Q  +  =  2%(0<K£  -  6n)  at  x  =  xu 


im (Oi>Q  -  =  0 


at  x  =  x2- 


(4.19a) 

(4.19b) 


To  recover  the  scattering  distribution,  we  must  solve  equation  (4.16)  with  the  mixed  boundary 
conditions  (4.19).  From  equations  (4.14)  and  (4.18),  it  follows  that 


r{9]p,6 in)  =  4>q(xi,p sinfl)  -  lg=ein  and  t(0;p,6in )  =  tPq(x2,P2(p)  sind) 


where  the  indicator  function  le=ein  equals  1  if  8  =  8m  and  it  equals  0  otherwise. 

If  the  potential  Vq(x, y)  is  constant  along  the  y-direction,  i.e.,  VQ(x,y)  =  Vq(x),  then  equa¬ 
tion  (4.16)  simplifies  to  the  separable  equation 

d2  - 

~  riiOi’Q  +  2 mVQ(x)i)Q  =  0.  (4.20) 

Since  Ex  =  p2(£t)/2m  is  simply  the  contribution  of  the  x-component  of  the  momentum  to  the 
kinetic  energy,  we  have  the  one-dimensional  Schrodinger  equation 


1  d2  *  * 

7T-  +  Vq(x)4>q  =  Ex^q 

2m.  oxz 


(4.21) 
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with  boundary  conditions  (4.19).  One  may  also  solve  the  boundary  value  problem  (4.21)  by  using 
the  transfer  matrix  method  presented  in  Section  3.1.1.  Note  that  since  the  solution  is  constant  in 
the  redirection,  the  semiclassical  impulse  force  is  normal  to  the  barrier  curve. 

To  solve  the  boundary  value  problem  (4.16)  for  general  potential  V ( x ,  y )  one  may  use  an  iterative 
solver 

Au(n+1)  =  Bu (n) 

where  utJ  is  the  discretization  of  V’cOeo  £j)>  A  is  the  finite  difference  operator 

Auij  =  -i(Ax)-2  [ui+1J  -  2 Uij  +  Ui-ij]  +  r)i(€j)uij  (4.22) 

and 

Biiij  =  —2  mJ:VijJ7~1Uij 

where  T  is  the  discrete  Fourier  transform  with  respect  to  y  and  Vj-j  =  V (. Xi,yj )  is  the  discretization 
of  the  potential. 

4.2.3  A  particle  method  for  the  semiclassical  Liouville  equation 

Following  initialization,  we  solve  the  Liouville  equation  using  the  particle  method  by  sampling 
N  particles  from  an  initial  distribution,  solving  Hamilton’s  equations  over  a  given  time  interval, 
and  then  fitting  the  data  to  an  appropriate  mesh.  By  linearity  of  the  Liouville  equation,  the 
particle  method  may  be  implemented  for  each  particle  independently,  permitting  us  to  speed  up 
computation  by  using  a  parallel  computer  cluster. 

A  particle  is  sampled  from  the  initial  distribution  deterministically,  for  which  we  associate  a 
weight 

wn  =  f  /(x,  p)  dx.  dp  (4.23) 

Jcn 

to  the  particle  (x°,p°,0),  or  using  Monte  Carlo  sampling,  for  which  we  associate  a  uniform  weight 
wn  =  A-1.  Monte  Carlo  sampling  is  important  in  higher  dimensions  because  it  mollifies  the 
“curse  of  dimensionality”  which  afflicts  deterministic  sampling,  restricting  it  to  a  rather  coarse 
mesh  in  higher  dimensions.  In  6-dinrensional  phase  space  (M3  x  M3),  a  billion  particle  sample  when 
sampled  deterministically  allows  only  30  grid  points  in  each  direction.  On  the  other  hand,  Monte 
Carlo  sampling  is  inefficient  for  nonstandard  distributions  and  the  solution  is  noisy  even  with  a 
substantial  sample  size. 
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At  each  time  step  we  use  a  symplectic  solver  [16]  to  compute  (xn+1,pn+  ,in+1)  where  we  take 
tn  =  nAt.  We  estimate  the  updated  position  of  the  particle 


x*  =  xn  +  Atpn  -  i(A l)2Vh(x")  (4.24a) 

p*  =  pn  -  i  A t  (W(xn)  +  W(x*))  .  (4.24b) 


If  x*  is  in  the  same  region  as  xn,  i.e.,  if  the  particle  has  not  crossed  the  barrier  Tg  during  the  time 
interval  [tn,tn+1],  we  let  (xn+1,p"+1)  =  (x*,p*).  If  x*  is  in  a  different  region  from  xn,  then  we 
determine  the  time  At*  of  the  barrier  crossing 


At* 


ci(xn) 

d(x*)  -  d(xn) 


At 


(4.25) 


where  d{x)  is  the  signed-distance  to  the  barrier.  The  time,  position  and  momentum  of  intersection 
with  the  barrier  are  estimated  by  the  solver  (4.24)  using  At* 

We  use  the  push  interface  condition  (4.2)  to  connect  the  appropriate  bicharacteristics  at  the 
barrier.  Since  the  bicharacteristics  are  not  unique,  either  a  Monte  Carlo  approach  or  a  deterministic 
branching  method  must  be  used  to  select  a  bicharacteristic  using  conditional  probabilities  based  on 
the  incident  momentum.  In  the  Monte  Carlo  method,  the  scattering  angle  is  chosen  by  randomly 
sampling  from  the  distribution  of  scattering  directions.  Once  an  outgoing  bicharacteristic  is  chosen, 
we  compute  the  position  (xn+1,pn+1)  of  the  particle  at  time  tn+1  by  using  the  solver  (4.24)  with 
the  remaining  time  step  given  by  At  —  At*  with  At*  defined  by  (4.25). 

The  unit  normal  vector  to  Tg  at  x*  may  be  calculated  either  analytically  or  approximated  by 
using 

n  =  Vd(x*)/|Vd(x*)| 


where  the  signed-distance  d(x*)  is  interpolated  linearly.  The  component  of  the  momentum  normal 
to  Tq  at  x*  is 

=  (p±,q±)  =  (p  •  n)n  =  {px\nx\  +Pj/|nJ,|)(|na:|,  |ny|). 

When  the  potential  V (x,  y )  is  continuous  along  the  length  of  the  barrier  curve,  there  are  only  two 
branches — one  transmitted  branch  and  one  reflected  branch.  In  this  case,  the  correction  to  the 
momentum  is 

(Ap,Aq)  =  ( p±,q± )  1  +  yjl  +  2?nAI4|p±|_2Nj  for  transmission, 

(Ap,  Aq)  =  — 2(|p_L|,  l^j) 


for  reflection. 
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Figure  4.1.  Phase  plane  bicharacteristics  and  the  associated  binary  tree.  By  considering  the  solution  in 
terms  of  a  binary  tree,  one  may  construct  a  deterministic  solution. 

In  the  case  of  two  branches,  it  is  convenient  to  use  a  deterministic  branching  algorithm  by  contin¬ 
uing  the  solution  along  both  transmitted  and  reflected  bicharacteristics.  Consider  a  binary  tree  with 
nodes  representing  the  intersections  of  the  bicharacteristic  with  the  barrier.  See  Figure  4.1.  To  each 
branch  we  associate  a  conditional  scattering  (transmission  and  reflection)  probability  Skj ■  We  track 
along  a  branch  using  the  solver  (4.24)  until  we  reach  a  new  node.  The  particle  information  (x,  p,  t ) 
is  saved  at  the  node  and  we  take  the  reflection  branch.  We  continue  in  such  a  manner  —taking  the 
reflection  branch  at  each  new  node — until  the  end  of  the  simulation  time.  The  probability  that  a 
particle  follows  the  kill  forward  global  bicharacteristic  is  the  product  of  the  conditional  probabilities 
Skj  for  each  node.  Therefore,  from  equation  (2.18)  the  contribution  is 

Nk 

wn,k  =  Wn  n  sk,j-  (4.26) 

3  =  1 

We  back  up  to  the  most  recent  node  that  has  an  unexplored  transmission  branch.  The  particle 
information  (x,  p,  t )  is  set  to  information  previously  stored  at  that  specific  node.  We  then  take  the 
transmission  branch,  continuing  as  above  until  the  end  of  the  simulation.  Once  all  transmission 
branches  have  already  been  explored,  i.e.,  once  we  have  backed  up  to  zeroth  node,  we  have  found 
all  the  forward  bicharacteristics  for  the  particle  initially  at  (x°,p°). 

The  solution  p(x,  t)  =  ff^°  /(x,  p,  t)  dp  is  reconstructed  by  interpolating  over  a  uniform  mesh 
using  a  smoothing  kernel  such  as  a  bicubic  spline.  Let  Ax  and  Ay  denote  the  mesh  spacing  and  let 
the  nearest  mesh  point  to  (. x,y )  be  ( Xi,yj )  for  some  Let  r  =  (x  —  Xi)/ Ax  and  s  =  (y  —  yj) / Ay 

denote  the  offset  from  that  meshpoint.  Since  we  are  interested  in  recovering  the  position  density, 
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we  do  not  need  to  reconstruct  over  the  momentum.  For  l,  rn  E  {—2, . . . ,  2}  define  mesh-constrained 
approximation  to  (x,  y )  as 

fi+ij+m  =  wn,kA(r,  l)A(s,  m)  (4.27) 

where 

A(r,  l)  =  a(r  +  l  +  |)  -  a{r  +  l  -  |) 

with  the  cut-off  function  [30] 

0  u  <  —  2 

^(2  +  u)4  — 2  <  u  <  — 1 

,  .  \  +  i(2u  —  rr3)  —  —  1  <  u  <  0 

|  +  |(2u  —  rt3)  +  |u4  0  <  u  <  1 

l-^(2  +  u)4  1  <  u  <  2 

1  2  <  u 

The  deterministic  and  Monte  Carlo  particle  methods  are  summarized  in  algorithms  4.2.3  and  4.2.3 
on  page  55. 

4.3  Numerical  Examples 

In  this  section,  we  present  two  examples  to  verify  the  numerical  scheme  and  validate  the  semi- 
classical  model.  Because  of  limitations  in  computer  resources  required  to  solve  the  von  Neumann 
equation,  even  using  an  indirect  method,  we  shall  limit  the  analysis  to  the  0(1)  Schrodinger 
wavepacket.  In  the  first  example,  we  consider  the  scattering  on  a  circular  step-potential.  This 
geometry  is  important  because  it  captures  phenomena  such  as  caustics  and  internal  reflection.  In 
the  second  example,  we  consider  the  scattering  on  an  electron  diffraction  grating  for  which  the 
potential  varies  on  the  quantum  length  scale  along  the  length  of  the  barrier.  Such  an  interface 
produces  multiple  scattering  angles. 

To  solve  the  time-dependent  Schrodinger  equation  we  use  a  pseudospectral  method  with  Strang 
splitting  similar  to  the  method  used  for  the  one-diminsional  Schrodinger  equation  in  Section  3.2. 
The  kinetic  and  potential  terms  are  split  so  that  for  each  time  step  we  have 

ip(x,  y,  t  +  At)  =  eAtB/2T~l  eAtAT  (eAtB/2^{x,  y,  tfj  (4.28) 

where 

A  =  exp(At— S—{kl  +  k2))  and  B  =  exp(Af— V(x,  y)) 

Omni  &  if? 
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Figure  4.2.  Reflection  (i?)  and  transmission  (T)  coefficients  (left)  for  the  absorbing  potential  V(x)  = 
?'|sech2(a;/4£)  (right).  Note  that  the  potential  is  “tuned”  to  absorb  the  momenta  near  p  =  1.06  used  in 
Example  4.3.1. 

and  the  operators  T  and  T~x  denote  the  two-dimensional  discrete  Fourier  transform  and  discrete 
inverse  Fourier  transform  with  respect  to  the  ( x,y )  and  (kx,ky)  variables.  When  the  potential  is 
discontinuous,  the  solution  exhibits  artificial  oscillations  unless  At  <  ( Ax)2 /e  and  At  <  e/V(x). 
Because  computer  memory  and  computing  time  constraints  in  two-dimensional  computation,  we 
weaken  the  conditions  that  we  used  in  the  one-dimensional  case  and  take  Ax  <  e/2,  allowing  us  to 
take  At  <  e/4. 

By  solving  the  Schrodinger  equation  over  a  periodic  domain  (rather  than  an  unbounded  do¬ 
main),  spurious  solutions  are  eventually  introduced  as  information  is  transmitted  across  the  domain 
boundaries.  By  embedding  the  domain  in  a  larger  domain  we  can  emulate  an  unbounded  domain 
for  a  sufficiently  short  simulation  time;  however,  this  approach  is  inefficient  especially  in  higher 
dimensions.  An  alternative  method  to  approximate  an  unbounded  domain  is  to  employ  an  absorb¬ 
ing  potential  VB(x,y )  near  the  boundaries  [18,  25].  By  adding  a  negative  imaginary  potential  that 
decays  rapidly  away  from  the  domain  boundaries,  we  have  the  modified  Schrodinger  equation 
d 

— y,  t )  =  \irn~1  Ai/j(x,  y,  t )  -  iV (x,  y)ip(x,  y,  t )  -  VB(x,  y)^(x,  y,  t) 

where  Vs(x,y)  >  0.  Such  a  potential  should  be  strong  enough  to  eliminate  (at  least  to  machine 
precision)  any  wave  information  passing  through  the  boundaries,  yet  not  so  strong  as  to  reflect  the 
wave.  In  addition,  the  potential  should  be  sufficiently  narrow  so  that  it  does  not  affect  the  solution 
away  from  the  boundary  or  require  an  overly  large  border.  In  this  case,  we  take 


VB(x )  =  Vosech2((x  —  Xb)/e£) 
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where  Xf,  is  the  position  of  the  domain  boundary,  Vo  is  the  barrier  strength,  and  l  is  the  characteristic 
barrier  width.  (One  may  examine  other  absorbing  potentials  by  using  transfer  matrices  discussed 
in  Section  3.1.1.)  For  this  potential  the  transmission  and  reflections  coefficients  may  be  found 
exactly  giving  [11] 


r(— ipl  —  7)r(— ipl  +  7  +  1) 
T(-ip£)T(-ipe  +  1) 


and  R 


T(—ipl  —  7)r(— ipi  +  7  +  1) 

r(— 7)r(7  +  1) 


where  7  =  —  ^  +  ^\/l  —  SiVol'2  and  p  is  the  normal  component  of  the  incident  momentum.  When 
pi  1,  7  ~  (1  —  i)V0  '  l,  and  hence  T  «  R  at  p  =  V0  .  By  adjusting  Vo  we  may  “tune”  the 
barrier  to  absorb  a  specific  range  of  energies.  Note  that  the  reflection  and  transmission  coefficients 
are  independent  of  e  and  by  taking  Ax  =  e/2,  we  may  specify  the  barrier  thickness  in  terms  of  grid 
points. 

To  compare  the  convergence  of  the  Schrodinger  to  the  semiclassical  limit  in  two  dimensions  we 
consider  the  following  LNerrors: 


the  L1-error  of  the  position  probability  density  function  (pdf) 

I P(x,y,t)  -  p(x,y,t)\  dydx 


the  L1  -error  of  the  cumulative  density  function  (cdf) 

rv 


p(r,  s,  t )  —  p(r,  s,  t)  ds  dr 


'  —OO  J  —OO 


dy  dx 


the  L1-error  of  the  marginal  probability  distribution  function  (mpdf) 


P(x,  y,  t )  -  p(x,  y,  t )  dy 


dx 


the  L1  -error  in  the  cumulative  marginal  distribution  function  (mcdf) 


/OO  px  /*oo 

/  /  p(s,y,t)  -  p(x,y,t) 

-00  J — 00  J — OO 


dy  ds 


dx. 


In  the  above  definitions,  p(x,  y,t)  =  f f^°  f(x,  y,p ,  q,  t)  dpdq  for  the  semiclassical  Liouville  solution 
and  p(x,y,t)  =  \ip(x,  y,  t)\2  for  the  Schrodinger  solution.  In  both  examples  we  take  the  effective 


mass  m  =  1. 
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4.3.1  Schrodinger  0(1)  wave  envelope  with  a  circular  barrier 


Consider  the  circular  barrier  with  unit  diameter 


V(x) 


0  x  6  fii 
\  x  6  Q2 


where  fli  =  {  x  |  |x|  >  1  }  and  1^2  =  •  Consider  initial  conditions 

n\  1  ~  x0)2  -  (y  -  y0)2\  fi(PoX  +  q0y)\  /Ann^ 

0)  =  715P exp  ( - - J“p( - 7 - )  <4'29) 

describing  a  symmetric  Gaussian  wavepacket  initially  located  at  (xq  .  yo )  traveling  with  momentum 
(po  ■  qo ) .  In  the  semiclassical  limit  we  take  the  initial  conditions 

f(x,y,p,q )  =  ^2  exp  — J°^2  — — 5(p-p0)8{q  -  q0).  (4.30) 


Let  (xo,yo)  =  ( — 1,  —  1) ,  (po)<Zo)  =  (§,  f )  and  a  =  j.  We  compute  over  a  square  domain  with  length 
L  =  4.  Using  equations  (3.3)  and  (3.7),  we  determine  the  reflection  coefficient  to  be 


Rip) 


for  a  particle  entering  fli  from  Q2 
for  a  particle  entering  O2  from  fli. 


Also,  using  equation  (4.5)  we  may  determine  the  critical  angle  6C  =  cos-1  (|p|-1)  •  For  the  initial 
speed  | po  |  ~  1.06,  the  critical  angle  0C  ~  0.34  and  the  momenta  of  the  particle  in  the  region  U2  is 
approximately  0.35. 

To  solve  the  Schrodinger  equation,  we  use  a  pseudospectral  method  with  Strang  splitting  (4.28) 
with  Ax  =  e/2  and  At  =  e/4.  To  mollify  spurious  reflections  and  transmissions  across  the  periodic 
boundary  conditions,  we  employ  an  absorbing  boundary  with  width  i  =  50e  =  lOOAx  encircling  the 
domain.  The  semiclassical  solution  is  computed  using  a  deterministic  particle  method  with  approx¬ 
imately  109  particles  and  reconstructed  using  a  mesh  spacing  Ax  =  100-1.  Since  the  semiclassical 
solution  is  reconstructed  over  a  coarser  mesh  than  the  Schrodinger  solution,  we  linear  interpolate 
the  semiclassical  solution  to  match  the  semiclassical  and  Schrodinger  solutions. 

The  marginal  probability  (position)  density  function  f  p(x,  y,  t)  dy  for  the  semiclassical  Liouville 
solution  and  the  Schrodinger  solution  for  several  values  of  e  are  shown  in  Figure  4.3  on  page  51.  Time 
evolution  of  the  Schrodinger  solutions  and  semiclassical  solutions  are  shown  in  Figures  4.4  and  4.5 
on  pages  56  and  57.  The  errors  in  the  two  solutions  are  listed  in  Table  4.1.  Based  on  our  study, 
we  find  the  convergence  rate  of  the  U-errors  to  be  about  first  order. 
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4.1.  Errors  in  solutions  of  Example  4.3.1  for  different  values  of  e 


£ 

50-1 

100-1 

200-1 

400^ 1 

convergence 

d-error  (pdf) 

4.52  x  1CT1 

2.90  x  10'1 

1.73  x  10_1 

1.24  x  10_1 

0.6 

l1  -error  (cdf) 

1.52  x  liT1 

6.01  x  10'2 

3.64  x  10”2 

1.72  x  10-2 

1.0 

d-error  (mpdf) 

3.20  x  10'1 

1.01  x  10'1 

5.03  x  10"2 

2.37  x  10-2 

1.2 

l1  -error  (mcdf) 

9.56  x  10~2 

3.08  x  10'2 

1.42  x  10”2 

7.20  x  10“3 

1.2 

Notable  phenomena  emergent  in  the  Schrodinger  solution  to  this  example  are  formation  of 
interior  caustics  and  internal  reflection.  See  Figures  4.4  and  4.5.  Suppose  a  particle  originally  in  a 
region  fli  with  potential  V\  is  transmitted  across  an  interface  Tq  near  the  critical  angle  and  enters 
a  convex  region  with  potential  Vi  >  V\.  The  particle  “creeps”  internally  along  the  interface 
Tq,  and  with  a  nonvanishing  probability  the  particle  is  trapped  in  the  region  of  higher  potential. 
While  the  semiclassical  model  accurately  captures  both  caustics  and  internal  reflection,  the  classical 
model  does  not. 


4.3.2  Electron  diffraction  grating 

Consider  the  semiclassical  potential 

v(x,y)  =  [V'i  if<^>ere 
1 0  otherwise 

where  Tq  is  the  y-axis.  Note  that  Tq  may  be  any  smooth  curve — we  take  the  Tq  to  be  y-axis  to 
simplify  analysis.  The  quantum  potential  Vq  =  lim£_o  Vq(£X'i  £y')  where  x'-axis  is  normal  to  Tq 
and  the  y'-axis  is  parallel  to  the  Tq.  Let  the  local  quantum  potential  (with  x  and  y  scaled  by  e) 
be  given  by 

Vq(x,v)  =  f(x)  (1  +  cosay)  with  f(x)  =  |(1  +  cos irx) 

if  x  £  [—  1, 1]  and  y£l  where  a  is  some  parameter.  Take  V{x,y)  =  0  elsewhere. 

We  begin  by  determining  the  scattering  coefficients  for  the  barrier.  Consider  a  particle  with 
momentum  p  (and  energy  E  =  \p2)  with  an  incident  angle  and  a  scattering  angle  0.  The 
y-component  of  the  incident  momentum  is  £;n  =  p  sin#in  and  the  y-component  of  the  scattered 
momentum  is  £  =  psinfk  Then  the  x-component  of  the  momentum  is  ??(£)  =  \Jp 2  —  £2.  Let  Vq 
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Figure  4.3.  Marginal  position  density  function  for  Example  4.3.1  for  e  =  (a)  50_1,  (b)  100-1,  (c)  200”1, 
(d)  400”1  at  time  t  =  2.  The  numerical  semiclassical  limit  is  indicated  by  *.  The  plots  are  offset  by  0.25  for 
clarity. 


be  the  Fourier  transform  of  'ijjQ  with  respect  to  y  as  defined  in  Section  4.2.2.  By  using  the  identity 
/  f(x)(l  +  cos  ay)ijj  Q(x,  y)e~liv  dy  =  ±f(x)  ('tpQ(x,^  +  a)  +  2  V>Q(a;,£)  +  ^q{x^  -  a) 
equation  (4.16)  becomes 

-q^2^q(x,0  -  V2(0'4’q(x,0  +  f(x)  (V’Q^^  +  a)  +  2$q(x,Q  +^q{x,£-  a) )  =0  (4.31) 


with  the  mixed  boundary  values 


d 


i'niOi’Q  +  =  2i??(0^(C  - 

/-  d  '• 

irjiOi’Q  ~  -q^Q-  =  0 


at  x  =  —  1 
at  x  =  +1. 


(4.32a) 

(4.32b) 


To  solve  the  boundary  value  problem  (4.31)  and  (4.32)  we  consider  a  finite  difference  method. 
Let  x-i  be  the  discretization  of  x  over  [—1,1]  using  m  grid  points  with  uniform  spacing  Ax.  Let 
be  the  discretization  of  £  using  a  uniform  spacing  =  a/ d  for  an  integer  d.  The  second-order 
centered-difference  approximation  of  (4.31)  is 


ZUij  +  Ui—\j 

(Ax)2 


VjUij  +  “1“  +  V'i,j—d)fi  —  0 


(4.33) 
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where  we  define  u.tj  =  'ipQ(xi,  7),  /)  =  f(xi),  and  ??j  =  r/(£j).  Note  the  convention  that  the  coeffi¬ 
cient  z  denotes  the  complex  constant  z  and  the  subscript  z  denotes  an  index.  As  a  simplification, 
one  may  restrict  £in  to  a  gridpoint  k  (£;n  =  7)  and  interpolate  over  the  scattering  coefficients  to  ap¬ 
proximate  intermediate  values.  The  second-order  approximation  of  the  boundary  conditions  (4.32) 
are 


U2,j  UO ,j  0 .  x 
lVjUlj  +  — 2Ax —  =  2ir]i6ik 


ITjjU. 


mj 


U'm  —  l,j 

2Ax 


=  0. 


(4.34a) 

(4.34b) 


Substituting  the  value  for  «0,j  from  the  left  boundary  condition  (4.34a)  and  substituting  the  value 
for  urn+i,-j  from  the  right  boundary  condition  (4.34b)  into  equation  (4.33),  we  have  the  equivalent 
system  of  equations 


((Ax)2??2  -  2(A xffi  -  2)  +  ui+1J  +  Ui-ij 

-  (Ax)2 fiUij+a  -  (Ax)2  =  Ofor  1  <  i  <  m 


(4.35a) 


2 U2,j  +  ((Ax)2??2  -  2(Ax)2  fi  -  2  +  i2(Ax)rjj)  u ij 

-  (Ax)2 fiUij+d  -  (Ax)2 fiuij-d  =  4i(Ax)r)j6jk 
2um-hj  +  ((A  xfrj2  -  2(A  x)2  fm  -  2  +  i2(Ax)rjj)  umj 


(4.35b) 


(4.35c) 


(Ax)  fm'U"m,j+d  (Ax)  fmUm,j—d  —  O' 

By  condition  (4.14)  we  have  that  Uij  =  0  if  |T;  |  >  p.  Furthermore,  the  solutions  Uij  =  0  if 
j  £  {. . . ,  k  —  d,  k,  k  +  d, . . .  }.  Therefore,  for  each  incident  momenta  we  solve  system  (4.35)  for 
j  e  {. . .  ,k  —  d,k,k  +  d, . . Let  { J}  be  the  ra-element  enumeration  of  this  set.  In  this  case,  we 
may  express  the  equations  as  the  system  Mv  =  b  where  the  nm-element  vector  v  is  defined  using 
v-i+mj  =  b  is  defined  using  bi+mj  =  4z(Ax)?/j,  and  M  is  the  block  tridiagonal  matrix 


(T<V  D 
D  T( 2) 


M  = 


\ 


\ 


D  T(n_1)  D  ! 
D  T(n7 


In  this  matrix,  D  are  m  x  m  diagonal  matrices  with  Dij  =  —  (Ax)2 fiSij  and  are  m  x  m 
tridiagonal  matrices 


T^p  =  ((Ax)2??2  -  2(A x)2U  -  2)  dij  +  8i+1>j  + 
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with  the  exceptions  T-p  =  (A x)2i)2j  —  2(A x)2fi  —  2  +  i2(Ax)r/j  for  i  =  l, m  and  T^J  =  2  and 


(J)  _ 


,(J) 

1  m.m—1 


From  (4.11),  the  transmission  coefficients  are  given  by  \vm+jn\2r]T  1r]j.  The  reflection  coefficients 
are  given  by  |1— t>i_|_jn|2  when  J  corresponds  to  k  incident  and  \vi+j\2rjpr]j  otherwise.  The  discrete 
scattering  angles  are  given  by 


9  j  =  —  sin 


— i  /sfc 

\P\ 


which  is  simply  the  Fraunhofer  diffraction  grating  formula 

J  A  =  d(sin#in  +  sin#/) 

with  wavelength  A  =  2ire/\p\  and  groove  spacing  d  =  2ira~1 . 

When  n  and  m  are  large,  we  can  use  a  few  manipulations  to  efficiently  solve  the  system  Mv  =  b. 
By  moving  the  top  row  of  M  to  the  bottom  (and  doing  the  same  for  the  vector  b ),  we  get  an 
equivalent  system  with  an  almost  Tipper  triangular  matrix 

/  D  T(2)  D  \ 


M'  = 


D  T^-1)  D 

£)  J'i.n) 

\TW  D  0  / 


To  invert  such  a  matrix  we  rely  on  the  useful  Sherman-Morrison- Woodbury  algorithm  [13] 

M’-1  =  {A  +  UVT)~1  =  A~l  -  A~lU{I  +  VtA-1U)-1VtA~1  (4.36) 


where  we  take 


(d  r(2)  d 


A  = 


\ 


D  r(n-T  D 
jj  y(«) 

I  ) 


U  = 


(°\ 


0 

W 


V  = 


fTA 

D 

0 

0 

v-v 


Because  A  is  upper  triangular,  it  is  easy  and  fast  to  invert  by  back  substitution.  The  matrix 
I  +  VT A~XJJ  is  a  dense  m  x  m  matrix  that  is  inverted  using  LJJ  decomposition.  The  number  of 
floating-point  operations  needed  to  invert  M'  when  n  and  m  are  large  is  10m2n+  | m3.  The  usual 
block  tridiagonal  algorithm  using  the  Thomas  algorithm  requires  ^m3n  floating-point  operations 
because  forward-substitution  does  not  preserve  the  sparsity  of  the  block  matrices. 
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If  the  incident  momentum  is  p  <  \a  then  n  =  1  and  for  each  incident  bicharacteristic  there 
is  one  transmitted  bicharacteristic  and  one  reflected  bicharacteristic.  In  general,  the  number  of 
transmitted/reflected  bicharacteristics  will  be 

N(£ in)  =  1  +  floor  [a~l\p  -  £in|]  +  floor  [a-1  \p  +  £in|]  <  2a_1p  +  1. 

To  validate  the  semiclassical  model,  we  took  a=  \  and  considered  the  initial  conditions  (4.29) 
and  (4.29)  with  a  =  1/16  and  (po,Qo)  =  (cos  (9in,  sin  0m),  (xo,yo)  =  0.3(—  cos  Q\n,  sin  9[n)  where 
6m  =  10°.  The  solve  the  semiclassical  model  can  be  solved  exactly  by  considering  the  method  of 
characteristics  and  the  scattering  coefficients  computed  numerically.  In  this  case 

P{x,y,t)  =  po{x  ,y  )  +  >  s(6'fc)po  [x - - — 

“  V  cos 

where  x*  =  x  —  tcos6[n  and  y*  =  y  +  t sin 9m  and  po(x,y)  is  the  position  density  of  the  initial 
distribution  (4.29). 

The  p(x,  y)  =  2  contours  of  the  position  density  for  the  Schrodinger  solution  and  the  semiclas¬ 
sical  solution  are  shown  in  Figures  4.6  and  4.7  on  pages  58  and  59  for  e  =  200-1  and  800-1.  The 
errors  in  the  two  solutions  are  listed  in  Table  4.2.  Although  the  solutions  have  roughly  first-order 
convergence  in  probability  density  functions,  the  solutions  fail  to  converge  in  the  cumulative  density 
functions.  Furthermore,  as  evident  in  Figures  4.6  and  4.7,  while  the  semiclassical  model  does  agrees 
with  the  Schrodinger  solution  for  small  scattering  angles,  there  is  a  significant  discrepancy  at  larger 
scattering  angles.  A  possible  explanation  for  this  discrepancy  is  underresolution  of  the  mesh  for 
the  Schrodinger  solution. 

4.2.  Errors  in  solutions  of  Example  4.3.2  for  different  values  of  s 


Oin  *  sin  +  sin  din  \  1 

x,y  cos  ek  x)  1(-cos^>°) 


£ 

1 

o 

o 

200-1 

400-1 

800-1 

convergence 

l1 -error  (pdf) 

6.09  x  10'1 

3.05  x  10'1 

2.25  x  10"1 

2.09  x  10”1 

0.8 

l1  -error  (cdf) 

4.86  x  10~2 

4.82  x  10~2 

5.00  x  10“2 

5.07  x  10"2 

— 

d-error  (mpdf) 

3.46  x  10'1 

1.81  x  10~2 

1.33  x  10_1 

1.04  x  10_1 

0.9 

l1  -error  (mcdf) 

6.57  x  10~2 

6.58  x  10'2 

6.82  x  10“2 

6.93  x  10“3 

— 
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Algorithm  4.1.  Deterministic  semiclassical  particle  method 

0.  Initialization.  Calculate  the  transmission/reflection  coefficients  associated  with  the  compo¬ 
nents  of  momentum  normal  to  the  interface.  Save  the  coefficients  in  a  table  over  which  to 
interpolate. 

1.  Choose  a  point  (x°,p0,t0  =  0)  and  calculate  the  weight  w  associated  with  the  initial  distri¬ 
bution  using  (4.23)  . 

2.  While  the  node  index  I  >  0 

(a)  Calculate  (x*,p*)  from  (xn,pn)  using  (4.24). 

(b)  If  x*  and  xn  are  both  in  the  same  regions,  take  (xn+1,pn+1)  =  (x*,p*).  Otherwise: 

i.  Compute  the  position  and  momentum  at  barrier  (x*,  p*)  using  (4.25).  Compute  the 
unit  normal  n  at  x*. 

ii.  Increment  the  node  index  I  and  save  (x*,p *,t*)  to  the  new  node. 

iii.  Take  the  reflection  branch  and  calculate  (xn+1,  pn+1)  using  (4.24)  with  timestep 
At  —  At*  given  by  (4.25). 

(c)  If  t  >  tmax 

i.  Reconstruct  the  solution  using  (4.27)  and  (4.26). 

ii.  Decrease  I  to  latest  node  with  an  unexplored  transmission  branch. 

iii.  Set  (x*,p*,t*)  to  value  stored  at  node  I. 

iv.  Take  the  transmission  branch  and  calculate  (xn+1,  pn+1)  using  (4.24)  with  At  — At*. 


Algorithm  4.2.  Monte  Carlo  semiclassical  particle  method 

0.  Initialization.  Calculate  the  scattering  distribution  associated  with  the  momentum  incident 
to  the  quantum  barrier.  Save  the  coefficients  in  a  table  over  which  to  interpolate. 

1.  Choose  an  initial  particle  (x°,p°)  from  the  initial  distribution  using  Monte  Carlo  sampling. 

2.  For  each  particle,  while  tn  <  tmax 

(a)  Calculate  (x*,p*)  from  (xn,p”)  using  (4.24). 

(b)  If  x*  and  xn  are  both  in  the  same  regions,  take  (xn+1,pn+1)  =  (x*,p*).  Otherwise: 

i.  Compute  the  position  and  momentum  at  barrier  (x*,p*)  using  (4.25)  and  compute 
the  unit  normal  ri,  at  x*. 

ii.  Use  Monte  Carlo  sampling  of  the  scattering  coefficient  s(9)  to  determine  the  scat¬ 
tering  momentum  p*. 

iii.  Calculate  (xn+1,pn+1)  using  (4.24)  with  timestep  At  —  At*  given  by  (4.25). 

3.  Reconstruct  the  solution  using  (4.27). 
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Figure  4.4.  Solutions  for  Example  4.3.1  for  e  =  50  1  (left)  and  e  =  100  1  (right)  at  times  t  =  0,  2,  4,  6,  8. 
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Chapter  5 

Extensions  to  the  Model  and  Further  Research 


In  this  chapter  we  consider  corrections  to  the  semiclassical  model  as  motivation  and  direction 
for  future  research.  Since  the  goal  of  this  chapter  is  to  explore  extrapolations  to  the  ideas  developed 
in  previous  chapters  as  a  vehicle  for  future  research,  the  results  are  presented  without  the  same 
level  of  rigor  as  in  previous  chapters. 

In  developing  the  semiclassical  model  we  imposed  several  limiting  assumptions  on  the  potential 
barrier.  Namely,  the  width  of  the  barrier  is  O(s),  the  distance  between  neighboring  barriers  is 
0(1),  the  change  in  the  potential  VxV(x)  is  0(1)  except  at  quantum  barriers,  and  the  coherence 
time  is  sufficiently  short.  Two  notable  potentials  that  are  exceptions  to  these  limitations  are  crys¬ 
talline  domains  and  mesoscopic  barriers.  Periodic  crystalline  domains,  such  as  the  Kronig-Penney 
model  [26],  consist  of  narrow,  closely-spaced  potential  barriers  or  wells.  Because  the  separation  be¬ 
tween  neighboring  barriers  is  O(e),  the  barriers  may  not  be  considered  independent.  Furthermore, 
since  the  barriers  extend  across  the  whole  the  domain,  the  classical  region  is  effectively  coupled  into 
the  quantum  region.  A  mesoscopic  barrier  arises  when  the  scaled  Planck  constant  s,  while  small, 
is  nonvanishing  or  when  the  barrier  width  is  substantially  larger  than  the  effective  support  of  the 
quantum  wavepacket.  In  this  case,  the  barrier  may  not  act  like  a  single  scat.terer  but  rather  like  a 
series  of  multiple  scatters. 

As  discussed  in  Section  2.3,  the  barrier  scattering  in  the  semiclassical  model  is  a  time- irreversible 
and  entropy-increasing  process  because  the  interface  condition  is  not  one-to-one.  That  is,  the 
interface  condition  combines  information  from  multiple  separate  bicharacteristics.  The  Schroding- 
er  equation,  on  the  other  hand,  is  time  reversible.  As  discussed  in  Section  3.1.1,  the  transmission 
and  reflection  coefficients  are  computed  by  solving  the  Schrodinger  equation  in  order  to  derive  a 
quantum  scattering  matrix  for  the  barrier 


S 


Q  = 


ri  t2 
t\  r2 


(5.1) 
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where  r3  and  tj  for  j  =  1,2  are  complex  numbers.  The  complex-valued  scattering  matrix  (3.5) 
contains  not  only  information  about  scattering  probabilities,  but  also  the  phase  shift  and  ultimately 
the  phase  delay  time.  The  semiclassical  scattering  matrix 

S=(n  S)  =  (glt.p  7$)  where  «.  =  vV  - ~2mVi 

discards  the  phase  information.  Since  we  only  track  the  particle  density,  the  solution  is  decoherent 
away  the  barrier. 

The  difference  between  the  semiclassical  solution  f(x,p,t )  and  the  Schrodinger  solution  ip(x,t) 
becomes  evident  by  examining  the  position  density  p(x,t).  The  Schrodinger  equation  and  the 
Liouville  equations  are  linear  with  respect  to  ip  and  /,  respectively.  But  the  position  density 
is  only  linear  with  respect  to  /.  Consider  ipi,  ip2,  fi  and  /2  with  p\  =  f  f\  dp  =  \ipi\2  and 
P2  =  f  f2  dp  =  \ip2 12-  The  position  density  for  the  superpositioned  solutions  f\  +  /2  is  p\  +  P2,  but 
the  position  density  for  the  superpositioned  solutions  ip\  +  ip 2  is 

Pi  +  P2  +  ^\/p\P2  cos(#i  -  02)  (5.2) 

where  6j  is  the  phase  of  ipj  for  j  =  1,  2.  The  probability  amplitude  ip\  +  ip2  contains  an  additional 
coherence  term  which  is  not  captured  by  the  semiclassical  interface  conditions  (2.16)  and  (2.17). 
Hence,  the  physical  observables  of  the  semiclassical  Liouville  solution  in  general  will  not  agree  with 
the  physical  observables  of  the  Schrodinger  solution.  This  point  was  highlighted  in  Section  2.3  for 
the  harmonic  oscillator  with  a  delta-function  barrier.  A  natural  correction  to  the  semiclassical  model 
is  one  that  not  only  tracks  density  along  the  bicharacteristics  but  also  phase-offset  information  from 
the  barrier. 

5.1  Coherent  semiclassical  model 

In  this  section  we  shall  consider  a  coherent  semiclassical  model  and  discuss  the  numerical  imple¬ 
mentation.  In  the  limit  as  the  scaled  Planck  constant  e  — >  0,  phase  offset  may  not  be  well-defined 
if  only  because  the  semiclassical  measurements  cannot  resolve  quantum  distances.  Furthermore, 
since  the  solution  is  less  likely  to  to  exhibit  coherence  in  two-dimensions,  we  limit  the  discussion 
to  one  dimension. 

Define  the  semiclassical  probability  amplitude  as 

$(x,P,t)  =  yj  f(x,p,t)el6{p) 
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where  0(p)  is  the  phase  offset  from  the  initial  conditions  &(x,p,  0)  =  \J f(x,p ,  0).  Define  the 
coherent  probability  density  as 

fcoh(x,p,t)  =  | $(x,p,t)\2  . 

If  $>(x,p,t)  satisfies  the  evolution  equation 

(f<f>  <9<3?  dx  dp 

dt  dt  dt  dx  dt  dp 

then  fcoh(x,p,  t)  satisfies 

<9/col  i  dx  <9/coh  dp  if. /'c0li  

dt  dt  dx  dt  dp 

Hence,  if  $(x,p,t)  is  a  solution  to  the  Liouville  equation  for  initial  condition  &(x,p,  0)  then 
fcoh(x,p,t)  is  a  solution  to  the  Liouville  equation  for  initial  condition  fCoh{x,p,0).  Furthermore, 
|3>i  +  d>2|2  =  fi  +  f 2  +  ZVfife  cos(6»i  —  62).  If  /coh  is  independent  of  p  or  if  /1  is  a  scalar  multiple 
of  f‘2 ,  then  this  condition  is  equivalent  to  condition  (5.2). 

When  there  are  several  bicharacteristics  stemming  from  a  scattering  barrier,  we  define  the 
semiclassical  probability  amplitude  as  the  superposition 

®(x,P,t)  =  ^2sk(H(x,p))$k{x,p,t) 

k 

where  &k(x,p,t)  is  the  solution  along  the  fe-th  bicharacteristic 

$k(x,P,t)  =  J  ${x,p,0)(pk(x,p,t;x,p)dxdp 
and  kth.  global  bicharacteristic  for  the  Hamiltonian  H(x,p)  is  defined  as 

c pk(x,p,t ;  x,p)  =  6(x(t)  -  x)6(p(t)  -p). 

We  define  scattering  term  Sk(H(x,p))  as  the  product  of  complex- valued  transmission  and  reflection 

coefficients  (5.1)  along  the  fcth  bicharacteristic.  The  definition  is  similar  to  the  superpositioning 

of  the  particle  density  of  the  semiclassical  model  Section  2.2,  but  in  this  case,  it  incorporates  the 

linearity  of  the  Schrodinger  solution.  Hence,  we  have  the  coherent  probability  density 

l  2 

fcoh(x,p,t)  =  \^2sk(H(x,p))$k(x,p,t)  . 
k 

We  define  the  position  density  in  the  usual  manner 


P(x,t) 


J  fcoh{x,p,  t)  dp. 
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Consider  implementation  of  the  coherent  semiclassical  model  numerically  using  either  the  de¬ 
terministic  or  the  Monte  Carlo  particle  method  discussed  in  the  preceding  chapter.  A  first-order 
finite  volume  approach  is  analogous  to  the  method  discussed  in  Chapter  3  may  also  be  constructed. 
In  the  deterministic  particle  approach,  we  take  the  weight 

Wk=  &(x,p,0)dxdp 

Jck 

for  a  cell  C At  each  node  j  record  either  Sj  =  t  or  Sj  =  r  where  t  is  the  comp  lex- valued 
transmission  coefficient  and  r  is  the  complex- valued  reflection  coefficient  defined  by  the  scattering 
matrix  (3.5)  for  the  barrier  associated  with  the  node.  Take  d>fc(xn,pn,  tn)  =  wk^(x  ~  xn)6(p  — 
pn)Y\skj.  Use  the  cutoff  function  (4.27)  over  phase  space  ( x,p )  to  associate  the  Dirac  measure 
5(x  —  xn)8(p—pn)  with  the  mesh  ( Xi,pj ).  Finally,  combine  the  results  with  the  previous  solution.  In 
the  Monte  Carlo  approach,  we  sample  N  particles  randomly  from  <h(x,p,  0)  and  take  wn  =  N~1//2. 
At  a  barrier,  define  the  transmission  and  reflection  probabilities  as  usual  T  =  \t\2  and  R  =  \r\2 . 
Sample  a  uniform  random  variable  X.  If  X  <  T,  the  particle  is  transmitted;  otherwise,  the 
particle  is  reflected.  Set  Wk  tWk/T  for  transmission  and  Wk  <—  rwk/R ■  for  reflection.  Take 
&k(,xn,pn,tn)  =  WkS(x  —  xn)5(p  —  pn).  Use  the  cutoff  function  (4.27)  over  phase  space  (x,p) 
to  associate  S(x  —  xn)8(p  —  pn)  with  the  mesh  ( Xi,pj ).  Combine  the  results  with  the  previous 
solution.  The  obvious  shortcoming  of  the  deterministic  approach  is  that  it  is  impractical  to  track  all 
bicharacteristic  branches  for  complicated  geometries  with  several  barriers,  such  as  in  a  crystalline 
material.  The  shortcoming  to  the  Monte  Carlo  approach  is  inherently  unstable  because  of  the 
rescaling. 

As  a  simple  example  we  compare  the  Schrodinger  equation,  the  semiclassical  model  and  the 
coherent  semiclassical  model  for  the  harmonic  oscillator  with  a  delta-function  barrier 

V  (x)  =  \x2  +  eaS(x) 


introduced  in  Section  2.3.  The  complex-valued  scattering  coefficients  for  the  delta-function  barrier 
are  [42] 

,  i'P  ,  a 

t  = -  and  r  = - . 

ip  —  a  ip  —  a 


64 


Let  a  =  \/3  and  take  the  Gaussian  initial  distribution 


f(x,P,  0) 
4>(x,p,  0) 


(7TCT2)-1/4  exp 
(7T6)_1  exp 
(7T6)-1/2  exp  [ 


(_x_xo)  exp 


(x  -  x0)2\ 

a2  y 


(p-Po)2\ 
e2^-2  ) 


(x  -  x0)2\ 
2a2  ) 


(p~Po)2\ 
2 s2a~2  ) 


with  e  =  0.01,  a  =  0.1,  xq  =  1  andpo  =  0.  The  semiclassical  model  is  solved  using  the  finite- volume 
method  developed  in  Section  3.1.2  and  the  coherent  semiclassical  model  is  solved  using  the  deter¬ 
ministic  method.  The  Schrodinger  equation  is  solved  using  the  Crank-Nicolson  scheme  (3.16)  with 
the  delta-potential  8{x)  approximated  by  (Ax)-1<5io-  The  coherent  semiclassical  model  accurately 
describes  the  coherent  behavior  exhibited  by  the  Schrodinger  solution.  The  solutions  are  plotted 
in  Figure  5.1  on  the  following  page. 


5.2  Conclusion  and  directions 


Chapter  2  mentioned  that  conservation  of  energy,  by  providing  a  constraint  on  the  solutions, 
creates  an  equivalence  class  of  bicharacteristics  but  does  not  yield  a  unique  solution.  An  interface 
condition  based  on  the  physical  characteristics  of  the  barrier  was  needed  in  order  to  isolate  appro¬ 
priate  mixing.  In  Chapter  3  and  Chapter  4  we  implemented  the  model  in  the  case  of  a  decoherent 
interface  condition.  For  the  class  of  one  and  two-dimensional  thin  barriers,  we  demonstrated  that 
the  Schrodinger  and  the  von  Neumann  solutions  converge  to  the  semiclassical  Liouville  solutions  at 
a  rate  of  0{e).  Finally,  we  introduced  a  coherent  semiclassical  by  modifying  the  interface  condition. 

While  the  coherent  semiclassical  model  is  valid  for  a  limited  set  of  problems,  the  existence  of  an¬ 
other  model,  consistent  with  the  Hamiltonian  preservation  principle  and  with  an  interface  condition 
that  is  as  physically  justifiable  as  the  decoherent  model,  underscores  the  importance  of  choosing 
the  appropriate  interface  condition.  A  wrong  interface  condition  can  only  give  a  wrong  solution. 
The  future  goal  of  this  research  is  to  overcome  the  shortcomings  of  the  coherent  semiclassical  model 
and  to  develop  a  robust,  viable  model  that  serves  as  the  logical  first  step  to  solving  problems  over 
periodic  and  mesoscopic  potentials. 


Schrodinger 


Decoherent  Semiclassical 


Coherent  Semiclassical 


Figure  5.1. 
the  example 


Comparison  of  Schrodinger,  decoherent  semiclassical,  and  coherent  semiclassical  solutions  to 
in  Section  5.1  at  t  =  0,  ir  ...  67t. 
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